AD-A228  490 


BJK  FILE  COPY 


Applied  Research  Laboratory 


Technical  Report 


OTIC 


ELECTE 
NOVI  4 1990 


pennState 


UlSTRIBUnON"STATEMENT  A 


Approved  'or  public  release; 


The  Pennsylvania  State  University 
APPLIED  RESEARCH  LABORATORY 
P.O.  Box  30 

State  College,  PA  16804 


THE  MODIFICATION  OF  BACKGROUND 
DENSITY  PROFILES  IN  THERMOHALINE  CONVECTION 

by 

J.  L.  Schramm 
H.  N.  Shirer 
R.  Wells 


Technical  Report  No.  TR  90-014 
November  1990 


Supported  by: 

Space  and  Naval  Warfare  Systems  Command 


L.R.  Hcttche,  Director 
Applied  Research  Laboratory 


Approved  for  public  release;  distribution  unlimited 


1 


\‘l  14® 


REPORT  DOCUMENTATION  PAGE 


form  «pprwd 
OAffl  No.  0704-0198 


1.  AGENCY  USE  ONLY  (laav«  bUnk) 


\  2.  REPORT  DATE 

November  1990 


4.  TITLE  AND  SUBTITLE 


S.  FUNOMG  NUMBERS 


THE  MODIFICATION  OF  BACKGROUND  DENSITY  PROFILES  IN 
THERMOHALINE  CONVECTION 

4-  AUTHORS 


J.  L.  Schramm,  H.  N.  Shirer,  R.  Wells 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESSEES) 

Applied  Research  Laboratory 

The  Pennsylvania  State  University 

P.  0.  Box  30 

State  College,  PA  16804 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

TR“  90-014 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESSEES) 

Space  and  Naval  Warfare  Systems  Command 

Department  of  the  Navy 

Washington,  DC  20363-5100 

10.  SPONSORING  /MONITORING 

AGENCY  REPORT  NUMBER 

N000039-88-C-005 1 

11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION/  AVAILABILITY  STATEMENT 

12b-  DISTRIBUTION  CODE 

Unlimited 

13.  ABSTRACT  (Minimum  200  wortni 


Density  profile  modification  in  thermohaline  convection  is  examined  with  a  sixteen- 
coefficient  spectral  model  derived  from  the  two-dimensional,  shallow  Boussinesq  equatioi s 
These  density  modifications  are  important  in  near-neutral  density  stratification^  becaus e 
they  can  significantly  alter  the  depth  and  location  of  the  near-surface  sound  channel  . 
in -the-- ©e-eem.  A  linear  stability  analysis  of  the  conductive  solution  gives  a  character:  s 
equation  from  which  the  Hopf  bifurcation  points  signalling  the  onset  of  convective, 
or  temporally  periodic,  solutions  can  be  determined.  Two  horizontal  wavenumbers  are 
included  to  allow  the  existence  of  a  double  Hopf  point,  which  when  separated,  may  yield 
diffeo-Hopf  points  on  the  periodic  solutions  and  so  allow  for  secondary  branching,  or 
torai,  solutions.  The  periodic  solutions  associated  with  the  Hopf  bifurcations  are 
the  torai  solutions  associated  with  the  diffeo-Hopf  points  are  then  studied. 


tic 


An  example  using  parameter  values  corresponding  to  those  of  melting  sea  ice  is  presentee 
In  this  case,  the  torai  solutions  are  always  unstable  owing  to  the  existence  of  a  diffec 
point  exclusively  on  the  unstable  periodic  solution  that  branches  second  from  the  conduc 
solution.  Thus,  the  competition  between  the  stabilizing  salinity  gradient  anH  HpgraMli 

14.  SURiECT  TERMS  (continued  on  next  page)  1 15.  NUMOER  OF  PAGES 


-Hop 

tive 

zing 


densitv  orofiles.  thermohaline  r.onver ti on  . 

Boussinesq  equation, 

Hopf  point,  melting 

sea  ice 

14.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

RWWSliified 

IB.  SECURITY  CLASSIFICATION 

ffiMed 

19.  SECURITY  CLASSIFICATION 

$?nrfl^n?ied 

20.  LIMITATION  OF  ABSTRACT 

SAR 

N5N  7540-01-180-5500 


Standard  form  198  <R*v  1-89) 

on  AHSt  Mo  ''<*  I 


THE  MODIFICATION  OF  BACKGROUND  DENSITY  PROFILES  IN  THERMOHALINE  CONVECTION 
Abstract  (continued) 


'temperature  gradient  allows  only  stable  steady  and  stable  periodic  solutions,- and- Hf 
CTs  theUt  solutions— ttraf  modify  the  background  density  profiles.  The  periodic  solutions 
are  characterized  by  an  unstable  density  profile  moving  alternately  between  the  interior 
and  the  boundaries.  Hence,  the  sound  channel  alternately  increases  and  decreases  in 
depth  and  in  location  relative  to  the  surface.  In  constrast,  the  steady  solutions  produce 
an  unstable  density  gradient  in  the  interior  of  the  domain.  Thus,  in  the  steady  case, 
there  is  always  a  near-surface  sound  channel. 


ABSTRACT 


iii 


Density  profile  modification  in  thermohaline  convection  is 
examined  with  a  sixteen-coefficient  spectral  model  derived  from  the 
two-dimensional,  shallow  Boussinesq  equations.  These  density 
modifications  are  important  in  near-neutral  density  stratifications 
because  they  can  significantly  alter  the  depth  and  location  of  the 
near-surface  sound  channel  in  the  ocean.  A  linear  stability  analysis 
of  the  conductive  solution  gives  a  characteristic  equation  from  which 
the  Hopf  bifurcation  points  signalling  the  onset  of  convective,  or 
temporally  periodic,  solutions  can  be  determined.  Two  horizontal 
wavenumbers  are  included  to  allow  the  existence  of  a  double  Hopf 
point,  which  when  separated,  may  yield  diffeo-Hopf  points  on  the 
periodic  solutions  and  so  allow  for  secondary  branching,  or  toral, 
solutions.  The  periodic  solutions  associated  with  the  Hopf 
bifurcations  and  the  toral  solutions  associated  with  the  diffeo-Hopf 
points  are  then  studied. 

An  example  using  parameter  values  corresponding  to  those  of 
melting  sea  ice  is  presented.  In  this  case,  the  toral  solutions  are 
always  unstable  owing  to  the  existence  of  a  diffeo-Hopf  point 
exclusively  on  the  unstable  periodic  solution  that  branches  second 
from  the  conductive  solution.  Thus,  the  competition  between  the 


stabilizing  salinity  gradient  and  destabilizing  temperature  gradient 
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allows  only  stable  steady  and  stable  periodic  solutions,  and  it  is 
these  solutions  that  modify  the  background  density  profiles.  The 
periodic  solutions  are  characterized  by  an  unstable  density  profile 
moving  alternately  between  the  interior  and  the  boundaries.  Hence, 
the  sound  channel  alternately  increases  and  decreases  in  depth  and  in 
location  relative  to  the  surface.  In  contrast,  the  steady  solutions 
produce  an  unstable  density  gradient  in  the  interior  of  the  domain. 
Thus,  in  the  steady  case,  there  is  always  a  near-surface  sound 
channel. 
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Chapter  1 
INTRODUCTION 

Acoustic  waves  constitute  one  important  mode  in  the  dynamics  of 
the  ocean.  Propagation  of  these  waves  is  greatly  affected  by  the 
salinity  and  temperature  fields;  the  vertical  profiles  of  salinity 
and  temperature  support  such  density-driven  processes  as  thermohaline 
convection,  salt-fingering,  and  small-scale  mixing.  Density 
differences  caused  by  these  processes  may  lead  to  the  development  of 
a  near-surface  sound  channel  that  traps  acoustical  energy  and  allows 
it  to  propagate  over  large  distances  (Apel,  1987).  By  understanding 
the  background  density  profiles  that  characterize  these  phenomena,  we 
can  gain  insight  into  how  acoustic  waves  propagate  in  the  vicinity  of 
these  density  variations. 

Focusing  on  thermohaline  convection,  we  consider  a  fluid  in 
which  both  the  temperature  and  salt  concentration  decrease  linearly 
with  height.  This  configuration  may  be  observed  in  estuaries  or  near 
melting  sea  ice.  In  thermohaline  convection,  the  two  components , 
heat  and  salinity,  act  in  opposite  ways  to  alter  the  vertical  density 
gradient.  Heating  from  below  or  cooling  from  above  acts  to 
destabilize  the  fluid  by  causing  the  density  of  the  fluid  at  the 
bottom  to  decrease  or  at  the  top  to  increase  to  values  sufficient  for 
convective  overturning  to  occur.  Freshening  above  or  salting  below 
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acts  to  increase  the  density  at  the  bottom  of  the  fluid,  thereby 
stabilizing  the  system. 

We  can  further  understand  the  physical  mechanisms  of 
thermohaline  convection  by  using  a  parcel  argument  (Huppert  and 
Moore,  1976).  If  a  parcel  is  initially  raised  from  its  equilibrium 
position  in  a  fluid  having  decreasing  temperature  and  salinity 
gradients,  then  it  comes  into  an  environment  that  is  cooler  and  less 
salty.  Because  temperature  diffuses  much  more  rapidly  than  salt,  the 
parcel  quickly  equilibrates  to  the  temperature  of  its  environment, 
making  the  parcel  saltier  and  more  dense  than  its  surroundings.  The 
parcel  sinks,  overshooting  its  initial  position.  Ihe  temperature 
again  equilibrates  rapidly  with  that  of  the  environment,  and  so  the 
parcel  becomes  fresher  than  its  surroundings.  It  is  therefore  less 
dense,  rises  toward  its  initial  position,  which  it  again  overshoots . 
An  oscillatory  response  occurs  that  can  be  either  (i)  damped,  or  (ii) 
periodic,  depending  on  how  fresh  the  parcel  is  as  compared  with  its 
environment.  If  the  gradients  driving  the  motion  are  small,  then  the 
parcel  is  only  slightly  fresher  than  its  environment,  and  so  it 
cannot  rise  to  its  initial  height.  Thus,  the  result  is  a  damped 
oscillation  to  its  equilibrium  position.  In  the  other  case,  the 
parcel  is  much  fresher  than  its  environment,  it  returns  to  its 
initial  position  and  a  stable  oscillation  occurs.  Based  on  the 
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conservation  of  mass,  vertical  motion  in  a  system  bounded  above  and 
below  will  lead  to  divergence  in  the  horizontal;  this  allows  the 
existence  of  two-dimensional  convective  cells. 

The  above  behavior  can  be  expressed  mathematically  in  the 
following  way.  For  a  given  salinity  gradient,  which  can  be 
represented  by  a  salinity  Rayleigh  number  Rag,  a  periodic  solution 

replaces  the  conductive  state  when  the  temperature  gradient, 
represented  by  a  thermal  Rayleigh  number  Raj,  exceeds  a  critical 

value  H^.  This  critical  value  is  known  as  a  Hopf  bifurcation  point 

from  which  emanates  a  periodic,  or  convective,  solution  characterized 
by  single  horizontal  and  vertical  wavenumbers .  This  periodic 
solution  is  only  one  of  the  infinite  nunber  of  branching  solutions 
represented  by  the  governing  equations  of  motion.  The  wavenumber  of 
the  periodic  solution  that  bifurcates  first  is  the  preferred  one  that 
is  expected  to  be  observed  in  the  fluid  (Laufersweiler,  1987). 

Veronis  (1965)  found  steady  and  periodic  solutions  in  thermohaline 
convection  using  a  model  consisting  of  five  coupled,  nonlinear 
ordinary  differential  equations.  In  the  Veronis  model,  the  value  of 
Raj  at  the  Hopf  bifurcation  point  is  much  less  than  the  value  of  Raj 

at  the  steady  bifurcation  point.  This  result  is  a  manifestation  of 
the  slow  diffusivity  of  salt.  In  this  model,  the  convective  flow  and 
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its  distortions  of  the  mean  temperature,  salinity,  and  velocity 
fields  are  minimally  represented. 

Using  the  Veronis  ( 1965 )  system  of  equations  with  moderate 
values  of  Rag,  Da  Costa  et  al.  (1981),  Knobloch  and  Proctor  (1981), 

and  Knobloch  et  al.  (1986)  found  oscillatory  solutions  ending  on 
unstable  steady  branches  as  the  thermal  forcing  rate  was  increased. 

For  large  values  of  Rag ,  they  found  a  period  doubling  cascade  to  a 

chaotic  solution  as  the  magnitude  of  Ra^,  was  increased.  This  is  one 

of  the  many  routes  to  chaos  that  have  been  proposed  (e.g.  Nese, 

1987). 

For  rotating  Rayleigh-Benard  convection,  Wells  and  Shirer  (1987) 
found  quasi -periodic ,  or  toral  solutions  that  involve  the  nonlinear 
coupling  of  two  separate  circulations  dominated  by  two  horizontal 
wavenumbers.  These  results  were  obtained  using  a  spectral  model  that 
contained  a  double  Hopf  point,  at  which  coincide  two  Hopf  bifurcation 
points  characterizing  periodic  solutions  of  two  different  horizontal 
wavenumbers.  Although  modification  of  the  background  state  was  not 
examined  for  the  toral  case,  the  model  of  Wells  and  Shirer  (1987)  did 
provide  an  explanation  of  cell  broadening  and  period  shortening  in 
rotating  convection.  Noting  that  thermohaline  convection  also  has 
double  Hopf  points,  we  investigate  whether  the  same  wavenumber 
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changing  mechanism  associated  with  toral  solutions  exists  in 
thennohaline  convection. 

Beginning  with  the  two-dimensional,  shallow  Boussinesq  system, 
we  convert  the  partial  differential  equations  to  a  dimensionless  form 
that  allows  us  to  define  Rag  and  Ratp  as  well  as  the  other 

dimensionless  parameters.  A  linear  stability  analysis  is  then 
performed  on  the  conductive  solution  to  these  equations  to  determine 
the  characteristic  equation  from  which  we  can  deduce  the  bifurcation 
points  signalling  the  occurrence  of  convective  solutions.  From  this 
analysis,  we  can  find  the  preferred  horizontal  wavenumber,  the  values 
of  Reup  for  the  steady  and  Hopf  bifurcations,  and  the  value  of  the 

aspect  ratio  a  and  the  magnitude  of  Rap  at  the  double  Hopf  point.  To 

find  approximations  of  the  solutions  branching  from  the  conductive 
solution  at  Hopf  bifurcation  points,  we  develop  a  nonlinear  spectral 
model. 

We  need  expansions  for  the  salinity,  temperature,  and  velocity 
fields  that  are  solutions  to  the  linearized  equations  of  motion  and 
that  satisfy  the  cyclic,  stress-free  boundary  conditions  for  flow  in 
a  rectangular  domain.  Here,  trigonometric  functions  suffice.  To 
capture  the  solution  at  the  double  Hopf  point,  we  must  include  at 
least  two  horizontal  wavenumbers,  m  and  j ,  and  vertical  wavenumber 


one  in  the  representations  of  the  dependent  variables.  We  may  think 
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of  a  double  Hopf  point  as  the  superposition  of  two  simple  Hopf 
points;  when  these  two  Hopf  points  are  separated,  a  toral  or 
secondary  branching  solution  may  occur  at  a  diffeo-Hopf  point.  Based 
on  the  results  of  Chang  and  Shirer  ( 1984 )  and  Wells  and  Shirer 
(1987),  we  find  that  representation  of  this  solution  requires 
catalytic  components  having  horizontal  wavenumbers  m+j  and  ro-j ,  and 
vertical  wavenumber  two.  These  catalytic  terms  represent  a  transfer 
of  stability  from  a  solution  having  one  horizontal  wavenumber  to  one 
having  another.  The  wavenumber  sets  (m,l),  (j,l),  (im+j,2),  and 
(m-j,2)  are  only  four  of  the  infinite  number  of  sets  involved  in 
solutions  to  the  partial  differential  equations. 

To  represent  the  changes  in  the  mean  background  salinity  and 
temperature  fields  by  the  convective  flow,  several  more  terms  in  the 
variable  expansions  are  needed  (Veronis,  1968).  The  interaction 
between  the  (m,l)  modes,  or  the  ( j ,  1 )  modes  can  be  represented  by 
components  having  horizontal  wavenumber  zero  and  vertical  wavenumber 
two.  To  account  for  the  modification  of  the  background  profiles  by 
toral  solutions,  we  add  temperature  and  salinity  coefficients  for 
horizontal  and  vertical  wavenumbers  zero  and  four,  respectively. 
These  (0,2)  and  (0,4)  terms  represent  the  attempt  by  the  circulation 
to  eliminate  an  unstable  density  gradient.  Summing  the  terms  in  the 


expansions  gives  a  total  of  sixteen  trigonometric  functions  needed  to 
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represent  toral  convective  solutions  and  the  modifications  of  the 
background  density  field  by  these  solutions. 

Using  this  sixteen-coefficient  low-order  spectral  model,  which 
is  composed  of  two  nonlinearly  linked  Veronis  models  representing 
periodic  solutions  having  horizontal  wavenumbers  m  and  j,  we  will 
examine  the  development  of  quasi -periodic  solutions.  We  study  this 
development  in  the  following  way.  Noting  that  the  occurrence  of  a 
diffeo-Hopf  point  is  related  to  the  separation  of  the  double  Hopf 
point  whose  value  depends  on  the  aspect  ratio  a,  we  choose  a  value  of 
a  near  that  of  the  double  Hopf  point  and  then  increase  the  value  of 
Raj.  past  that  for  the  first  Hopf  bifurcation  point.  The  resulting 

stable,  periodic  solution  is  characterized  by  either  wavenumber  m  or 
j,  and  we  increase  the  value  of  Raj  further  until  a  diffeo-Hopf 

bifurcation  occurs  on  this  periodic  solution.  This  diffeo-Hopf 
bifurcation  point  can  be  found  using  a  stability  analysis  of  the 
periodic  solution.  At  the  diffeo-Hopf  bifurcation,  temporal  behavior 
is  excited  in  the  catalytic  terms,  and  so  a  quasi -periodic  or  toral 
solution  becomes  possible.  We  are  interested  in  the  thermohaline 
circulations  that  correspond  to  these  toral  solutions  and  also  in  the 
way  these  quasi -periodic  solutions  modify  the  background  density 
profiles  (Veronis,  1968). 

We  begin  in  the  following  chapter  by  developing  the  model.  In 
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Chapter  3,  we  construct  a  system  of  perturbation  equations  that  will 
be  used  to  determine  the  stability  of  the  periodic  solutions.  This 
examination  of  stability  will  allow  us  to  locate,  via  an  eigenvalue 
analysis,  the  diffeo-Hopf  bifurcations,  as  well  as  other  ones  — 
harmonic  bifurcations,  and  period  doubling  bifurcations  —  on  the 
periodic  solutions.  The  results  of  this  analysis  and  some  examples 
of  the  circulations  corresponding  to  steady  and  periodic  solutions 
are  examined  in  Chapter  4.  In  these  examples,  parameter  values 
corresponding  to  those  characterizing  various  oceanic  flows  are  used. 
Discussion  focusing  on  the  modification  of  the  background  density  by 
these  steady  and  periodic  solutions  is  also  included.  Concluding 


remarks  are  found  in  Chapter  5. 
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Chapter  2 
MODEL  DEVELOPMENT 

In  double-diffusive  convection,  several  types  of  motion  are 
possible.  The  flow  resulting  from  a  particular  mechanism  depends  on 
whether  the  driving  energy  comes  from  the  substance  of  higher  or 
lower  diffusivity  (Turner,  1973).  In  thermohaline  convection,  the 
diffusivity  of  temperature  is  100  times  that  of  salt.  This  large 
difference  in  diffusivities  produces  a  destabilizing  process  that 
maintains  convective  motion  in  both  the  diffusive  regime  and  the  salt 
finger  regime. 

Here,  we  focus  on  the  diffusive  regime  *r.  which  cold,  fresh 
water  overlies  warm,  salty  water,  inis  situation  can  be  produced  in 
a  layer  of  fluid,  for  example,  by  heating  aijd  salting  it  from  below. 
In  this  case,  a  stabilizing  salinity  gradient  and  a  destabilizing 
temperature  gradient  are  created. 

The  diffusive  mode  may  occur  at  the  mouths  of  rivers  or  near 
melting  ice  islands  in  the  Arctic  Ocean  (Turner,  1973).  Indeed,  sea 
ice  strongly  influences  che  suppression  or  development  of  convection. 
For  example,  melting  sea  ice  increases  the  static  stability  of  a 
column  by  adding  fresh  water  to  the  surface,  while  the  formation  of 
sea  ice  destabilizes  the  water  column  by  introducing  salt  into  the 


upper  layer. 
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To  study  the  circulation  induced  in  the  diffusive  regime  by  the 
temperature  and  salinity  gradients,  we  introduce  the  Boussinesq 
equations  of  motion.  A  linear  stability  analysis  of  the  conductive 
state  reveals  some  properties  of  flow  in  the  diffusive  regime. 
Finally,  a  Fourier  series  solution  to  these  partial  differential 
equations  is  assumed,  and  then  truncated  to  produce  a  low-order 
spectral  model,  as  shown  in  the  final  section. 

2.1  The  Boussinesq  Equations  and  Boundary  Conditions 

In  the  Boussinesq  system,  variables  can  be  expressed  as 
sums  of  the  conductive  and  convective  elements: 


y(x,z,t) 

=  V*  (x,2 

:,t)  , 

(2.1) 

PlX,Z,t) 

^00  * 

pQl z)  +  p’ (x,z,t)  , 

(2.2) 

T(x,z,t) 

=  Vz> 

+  T’ (x,z,t)  , 

(2.3) 

S(x,z,t) 

=  Vz> 

+  S’ (x,z,t)  , 

(2.4) 

p(x,z,t) 

=  p0(z) 

+  p' (x,z,t)  , 

(2.5) 

where  v,  p,  T,  S,  and  p  are  velocity,  density,  temperature ,  salinity, 

and  pressure.  The  subscripts  o  and  oo  denote  the  conductive  state 
and  the  primes  denote  the  convective  state.  Here,  we  note  that  a 


three-variable  density  decomposition  is  used  in  (2.2)  instead  of  the 
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usual  two-variable  one.  Substituting  (2.1)- (2. 5)  into  the  full 
governing  thermohaline  equations  gives  the  motionless  and  conductive 
state.  This  conductive,  motionless  state  is  a  hydrostatic, 
horizontal  layer  of  fluid  of  depth  and  density  pQQ  +  PQ(z)  and  is 

characterized  by  a  quadratic  profile  of  pressure  and  by  linear 
vertical  profiles  of  density  pq,  temperature  TQ ,  and  salinity  So; 

these  are  given  by 

3p0/3z  =  -g  <p0  +  p0fl)  ,  (2.6) 

To  =  Too  "  z  AT/^T  -  t2-7> 

S0  =  SQ0  -  z  AS/zt  ,  (2.8) 

where  poq,  pQ0,  Tqo,  and  SQ0  are,  respectively,  the  magnitudes  of 

density,  pressure,  temperature,  and  salinity  at  z  =  0.  The  vertical 
differences  of  S  and  T  are  given  by  AT  and  AS.  In  the  diffusive 
regime,  the  values  of  AT  and  AS  are  both  positive.  The  relationship 
between  pQ( z),  TQ(z),  and  SQ(z)  is  given  below. 

Restricting  motion  to  the  two  dimensions,  x  and  z,  we  may  follow 
Veronis  ( 1965 )  and  write  the  Boussinesq  equations  as 

9v’ 

+  v* •  Vv*  ~  Vp’~  gk(p  +  p’)/pqo  +  vV*v’  ,  (2.9) 

~  '  *00  ~  ~ 


V  •  V*  =  0  . 


(2.10) 
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The  conservation  equations  for  temperature  and  salinity  are 

|1’  +  v>  •  VT’  -  =  Kp^T'  ,  (2.11) 

—  1 

H’  +  v’  •  VS'  -  =  KgV*S ’  .  (2.12) 

~  "T 

Here,  the  thermal  and  solutal  diffusivities  are  given  by  Kp  and  Kg, 

and  v  is  the  kinematic  viscosity. 

Since  the  density  expressed  in  (2.2)  and  (2.9)  depends  on  both 
temperature  and  salinity,  we  can  write 


p 


oo 


(2.13) 


as 

P  =  P00  (1  -  «To  +  PSQ  -  aT’  +  PS’)  ,  (2.14) 

where  the  coefficients  of  expansion  are  given  by  «  =  1/T00  and  P  = 
1/SQ0  (Turner,  1973).  Combining  (2.2)  and  (2.14)  and  equating  the 
background  and  convective  parts,  we  get 


P0  =  P00  <-«T0  +  PS0)  =  P00(«hT  -  0hS)  z/zp  ,  (2.15) 

end 

p'  -  pQo  ( -<xT’  +  PS’)  ,  (2.16) 
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where  the  linear  relation  for  pQ  is  found  by  substituting  (2.7)  and 
(2.8)  for  SQ  and  TQ. 

By  introducing  the  stream  function  V,  where  u’  =  -3V7&Z  and  w’  = 
'di^/'dx,  and  cross -differentiating  (2.9)  to  eliminate  pressure,  we 
obtain  the  vorticity  equation 

I^V2*)  +  Jbfr.V2*)  =  g[«  |J’-  P  If’]  +  .  (2.17) 

The  thermodynamic  equations  for  heat  and  salt  can  be  rewritten  as 

W*  =  Kr^T’  *  (2-18) 

J(*,S’>  -  JS  ft  =  Kg^g,  f  ( 2 . 1.9 ) 

df  ag  af  ag 

where  the  Jacobian  is  defined  as  J(f,g)  =  dz  ~  dz  dx’  with  this 
notation,  we  can  write  the  boundary  conditions  as 

V  =  =  T’  =  S’  =  0  at  z  =  0,  Zj  .  (2.20) 

These  constraints  imply  stress-free  boundaries  that  are  perfect 
conductors  of  heat  and  salt.  In  the  horizontal,  cyclic  boundary 
conditions  ere  assumed  at  x  =  0,  L. 
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2.2  Dimensionless  Forms 


We  assume  that  the  horizontal  domain  is  infinite,  and  that  the 

characteristic  length  in  the  x-direction  is  L.  To  simplify  the 

Boussinesq  system,  we  introduce  the  following  dimensionless  forms: 

v  -  *  %  , 

(2.21) 

x  =  x*L/'2tt  , 

(2.22) 

*  . 

Z  =  Z  Z , 

(2.23) 

t  =  t*Zj  L/2nz  fcj.  , 

(2.24) 

T’  =  T*v  ^  T00  *3/g  , 

(2.25) 

S’  =  S*I/  Ky  S00  n3/g  Z^3  , 

(2.26) 

where  the  starred  variables  are  dimensionless.  An  aspect  ratio  a  may 


be  defined  as 


a  =  2  &j,/L  ,  (2.27) 

a  Frandtl  number  as 

P  =  v/Krp  ,  (2.28) 

and  a  ratio  r  of  the  diffusivities  as 


t  ~  Kg/Kj  . 


(2.29) 
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To  represent  the  forcing  applied  to  the  system,  we  define  a 


salinity  Rayleigh  number 


Rag  = 


*  «T  Soo  ** 


(2.30) 


and  a  thermal  Rayleigh  nunber 


Raj.  = 


^  %  Too  n 


(2.31) 


Positive  values  of  Raj  and  Rag  represent  a  destabilizing  temperature 

gradient  and  a  stabilizing  salinity  gradient,  respectively,  as  found 
in  the  diffusive  regime.  Here,  both  Rag  and  Raj  are  defined  using 

(Oj.  This  will  present  no  difficulties  in  interpreting  the  results, 

and  will  also  make  the  pQ  -  constant  line  correspond  to  the  Raj  =  Rag 

line  (Section  2.3).  The  dimensionless  Laplacian  operator  is  given  by 


^<  ,  =  a2  +  alL-t  . 


3x*2  3z*2 


Writing  (2.17)— (2.19)  in  dimensionless  form  gives 


(2.32) 


-j(vV)  +  J(y*,W*>  =  Pf^  -  |  vV  , 

c>  [_9x  3x  J 


(2.33) 


^  +  J(y* ,T* )  -  Raj  ^ 
dt  1  dx 


br*  _  V2!* 


(2.34) 


^  +  J(V*,S*>  -  Rag  ^  =  a  ^s* 

at*  ^  ax  a 


(2.35) 
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The  above  introduction  of  dimensionless  variables  enables  us  to 
specify  completely  the  system  by  knowing  the  values  of  the  five 
dimensionless  parameters  P,  a,  Rap,  Rag,  and  r.  For  the  diffusive 

regime,  we  study  the  case  in  which  the  values  of  both  Rayleigh 
numbers  are  positive.  Here,  P  and  r  characterize  the  properties  of 
the  fluid,  while  Rap  and  Rag  represent  externally  applied  forcing 

rates  (Huppert  and  Moore,  1976).  From  the  properties  of  salt  water, 
we  find  that  P  =  7  and  r  =  0.01.  Thus  the  problem  is  reduced  to  one 
in  the  three  parameters  Rag,  Rap,  and  a. 

2.3  Linear  Stability  Analysis 

Expressions  characterizing  some  of  the  properties  of  the 
diffusive  regime  can  be  found  using  a  linear  stability  analysis  of 
the  conductive  state.  Because  this  state  is  motionless,  the  Jacobian 
terms  from  ( 2 . 33 )-( 2 . 35 )  are  dropped,  producing  the  linear  equations 


dt 


*<vV)  =  Pp^f  -  asyl  +  f  W* 
*  [dx  dx  J  a 


at  1  ax  a 


K  -  **  K  =  I 

at  ®  ax  a 


(2.36) 

(2.37) 


(2.38) 


We  assure  solutions  of  the  form 


>*  ~  eXtsin(kx)sin(nz)  , 


T*  ~  eXtcos(kx)sin(nz)  , 


S*  ~  eXtcos(kx)sin(nz)  , 


(2.39) 


(2.40) 


(2.41) 


where  k  is  a  horizontal,  integral  wavenumber,  and  n  is  a  vertical, 
integral  wavenumber.  Also,  the  eigenvalue  X  is  a  complex  number 
whose  real  part  represents  the  growth  rate  and  imaginary  part 
repesents  the  oscillatory  behavior.  Subs i tut ing  (2.39)-(2.41 )  into 
( 2 . 36 )  —  ( 2 -  38 )  leads  to  the  following  characteristic  equation: 


x3+  ^-(l+P+T)  -  j^Pka(Raj.-Rag)  -  ^(P+tP+t)]  +  ^(-Ra^+Rag )  + 

=  0  ,  (2.42) 

a3 

where  <  =  a2k2  +  n2.  When  X  satisfies  this  cubic  equation,  a 
solution  to  ( 2 . 36 )  —  < 2 . 38 )  of  the  form  (2.39)-(2.41 )  exists.  Of  the 
three  roots  to  (2.42),  one  will  always  be  real .  The  other  two  roots 
may  be  either  real  or  complex  conjugates  (Baines  and  Gill,  1969). 

To  find  the  bifurcation  point  for  the  steady  solution,  where 
stability  is  lost  by  a  negative  eigenvalue  becoming  positive,  and 


therefore  going  through  zero,  we  substitute  X  =  0  into  (2.42)  to  give 
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Here,  Ra^,c  represents  the  onset  of  an  infinite  number  of  steady 

solutions,  one  for  each  combination  of  horizontal  and  vertical 
wavenumbers .  Exchanges  of  stability  usually  occur  at  such  critical 
values  of  the  forcing  parameter,  and  these  critical  values  sue 
determined  by  other  parameters  of  the  system.  We  can  minimize  the 
value  of  Ra^c  to  find  which  parameter  values  give  the  smallest 

forcing  rate  required  to  balance  energy  input  and  dissipation 
( Laufersweiler ,  1987 ) . 

The  existence  of  a  minimum  value  for  Rapc  indicates  that  a 

preferred  response  is  possible  in  the  fluid.  At  a  value  of  Ra^  less 

than  its  minimum  critical  value,  perturbations  of  all  wavelengths 
will  damp  because  every  eigenvalue  is  negative.  At  a  value  of  Rap 

greater  than  the  minimum  critical  value,  some  perturbations  are  able 
to  grow.  The  perturbations  that  first  grow  give  the  preferred  mode 
and  are  associated  with  a  certain  aspect  ratio. 

To  minimize  the  magnitude  of  Rapc,  we  differentiate  (2.43)  with 

respect  to  a2k*  to  get 

RaItain  =  “7s  +  H  *  <2,44) 

Here,  Rapmn  is  the  minimum  value  of  Rapc  for  the  development  of 

steady  convective  flow  and  occurs  when  a2k2  =  1/2  and  n  =  1 . 
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Physically,  this  corresponds  to  the  minimum  vertical  temperature 
gradient  required  to  support  convection.  We  note  that  setting  Rag  = 

0  in  (2.43)  and  (2.44)  gives  the  same  expression  for  Ra«|,c  and  Ra^^^ 

as  those  in  Rayleigh-Benard  convection. 

Since  r  <<  1,  the  first  term  in  (2.43)  dominates  when  Rag  > 

0.0675,  or  for  almost  all  values  of  Rag*  This  observation  implies 

that  the  destabilizing  temperature  gradient  must  exceed  the 
stabilizing  salinity  gradient  by  a  factor  of  100  before  steady 
convection  can  occur  . veronis,  1965).  However,  satisfying  this 
requirement  v  *i'  ;  lead  to  a  vertical  density  profile  that  is  very 
unstable  gravitationally.  A  bifurcation  from  a  trivial  solution 
directly  to  a  steady  one  would  require  too  large  a  difference  between 
the  temperature  and  salinity  gradients.  Intuitively,  this  result 
seems  aphysical,  leading  us  to  ask  whether  other  bifurcations  are 
possible. 

To  find  the  transition  from  the  conductive  solution  to  a 
convective  solution  via  a  Hopf  bifurcation,  we  let  X  r  ixm,  where  Am 

is  real.  Substituting  this  into  (2.42)  gives  a  complex  cubic 
polynomial  having  real  part 

E^(Ras-rRaT)  +  =  0  , 

a 


"xm  4-(1+P+t>  + 


(2.45) 
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and  imaginary  part 

“Xmf  ‘  Xm[Hc2(RaT"RaS)  "  ^(p+Tp+T)]  =  0  •  (2.46) 

To  obtain  the  Hopf  bifurcation  points,  we  solve  (2.46)  for  Am2  and 

note  that  *ffi2  >  0  must  hold  to  guarantee  the  existence  of  the  Hopf 

point.  Putting  (2.46)  into  (2.45)  gives  the  following  expression  for 
H^,  the  value  of  Ra^,  at  the  Hopf  bifurcation  point: 


Rag(P+T)  ( )  3  (  1+t+t/P)  (  1+t+P)  -  r 

1+P  +  a*k*(l+P) 


(2.47) 


Here,  (2.47)  represents  the  onset  nf  an  infinite  nunber  of  periodic 
solutions,  one  for  each  combination  of  horizontal  and  vertical 
wavenumbers.  Since  r  =  0.01,  and  P  =  7,  (2.47)  can  be  well 
approximated  by 


■V  =  nfr  ♦  ■  <*•«»' 

To  find  the  Hopf  bifurcation  point  1^  at  which  the  conductive 

regime  undergoes  a  transition  to  a  periodic  regime  of  horizontal 
wavenumber  m  and  vertical  wavenumber  1,  we  set  k  =  m  and  n  =  1  in 
(2.47).  For  the  Hopf  point  H-  denoting  the  onset  of  a  periodic 

solution  of  horizontal  wavenumber  j  and  vertical  wavenumber  1 ,  we  set 
k  =  j,  and  n  =  1.  Figure  1  shows  how  the  values  of  HtQ  and  Hj  vary 
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Figure  1  The  dependence  of  the  Hopf  bifurcation  points  and 

the  aspect  ratio  where  Rag=65,  m=2,  j=l ,  P=7,  and  r=0.01 
These  two  curves  intersect  at  the  double  Hopf  point 
given  approximately  by  a^rO. 49343. 


22 


with  a  for  the  case  of  m  =  2  and  j  =  1 .  In  a  system  having  two 
horizontal  wavenumbers  m  and  j,  the  double  Hopf  point  occurs  when 

the  Hopf  points  and  Hj  coincide  and  m-periodic  and  j-periodic 

solutions  are  created  simultaneously.  The  aspect  ratio  a^  at  which 

the  double  Hopf  point  occurs  is  given  by  setting  k  =  m  and  k  =  j  in 
(2.47)  and  then  equating  the  resulting  right  sides  to  give 


Vi 


2  _ 


m2/3  -  j2/a 
m2j2/3  -  j2m2/3 


(2.49) 


We  note  that  this  expression  for  a^  matches  that  for  double  steady 


bifurcation  in  Rayleigh-Benard  convection  (Hirschberg,  1987).  For  m 


=  2  and  j  =  1 ,  we  find  that  a^  =  0.49343.  This  point  is  of  special 

interest  since,  if  a  toral  solution  is  present,  then  it  appears  when 
the  double  Hopf  point  is  separated  into  two  single  Hopf  points  by 
increasing  or  decreasing  from  a^j  the  value  of  the  aspect  ratio  a 

(Figure  1).  We  note  that  the  five-coefficient  Veronis  (1965)  model 
Vm  or  V j  of  wave number  m  or  j  is  not  capable  of  supporting  a  double 

Hopf  bifurcation,  or  a  toral  solution,  since  a  toral  solution 
requires  two  horizontal  wavenumbers,  and  Vm  or  Vj  contain  only  one. 

In  the  absence  of  diffusion  and  viscosity,  there  is  no  Hopf 
bifurcation,  and  the  vertical  density  distribution  alone  dr  ermines 
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the  response  of  the  fluid  to  small  displacements  (Figure  2).  Above 
the  line  Rag  =  Ra^,  denoting  the  isopycnal  case,  a  parcel  displaced 

upwards  will  be  lighter  than  its  surroundings  and  its  motion  will 
exhibit  exponential  growth  away  from  its  initial  position.  Below 
this  line,  where  density  decreases  with  height,  parcels  will 

»»  It  II 

oscillate  at  the  Brunt-Vaisala  frequency.  the  constant  density 
line,  a  displaced  parcel  will  have  the  same  density  as  its 
surroundings  and  so  will  remain  in  its  new  position. 

Adding  diffusion  and  viscosity  produces  Hopf  bifurcation  points 
(2.47)  (Figure  3).  This  Hopf  bifurcation  line  is  an  adjustment  to 
the  constant  density  line.  Now  there  is  growth  to  a  periodic 
solution  above  the  Hopf  line,  and  damped  oscillations  below  this 
line.  Adding  dissipative  forces  allows  damped  oscillations  when 
density  increases  with  height  at  parameter  values  where  exponential 
growth  occurred  in  the  inviscid  case.  In  contrast,  there  are  still 
periodic  solutions  at  parameter  values  where,  in  the  inviscid  case, 

there  were  oscillations  at  the  Brunt-Vaisala  frequency. 

The  Hopf  and  steadv  bifurcation  boundaries,  as  well  as  the  pQ  - 

constant  line  given  by  Ra^Rag,  are  shown  in  Figure  4.  In  regions  A 

and  B  (shaded)  of  this  figure,  there  is  no  development  of  motion,  and 


ail  perturbations  damp  to  the  stable  conductive  solution.  In  region 


Figure  2  Convective  regimes  in  the  inviscid  case  where  positive 
values  of  Ra^  and  Rag  represent  destabilizing  and 

stabilizing  gradients,  respectively.  Displaced  parcels  i 
the  regime  above  the  isopycnal  line  Rag=Rap  will  grow 

exponentially  away  from  their  initial  positions,  while 
those  in  the  regime  below  the  Rag=RaT  line  will  oscillate 

««  •«  ft 

at  the  Brunt- Vaisala  frequency. 
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Figure  3  Convective  (unshaded)  and  conductive  (shaded)  regimes  when 
there  are  diffusion  and  viscosity.  Adding  these 
dissipative  forces  creates  Hopf  bifurcation  points,  which 
now  separate  the  regimes .  In  the  inviscid  case  ( Figure  2 ) , 
the  constant  density  line  divides  the  two  regimes.  The 
shaded  area  represents  damped  oscillations  toward  the 
conductive  solution.  Note  that  in  some  cases,  adding  these 
forces  allows  damped  oscillations  to  occur  even  when 
density  increases  with  height. 


Figure  4  Schematic  depiction  of  the  convective  (periodic),  salt 

finger,  and  conductive  regimes.  The  shaded  area  represents 
the  stable  conductive  solution.  The  Hopf  and  steady  lines 
denote  periodic  and  steady  bifurcation  points, 
respectively,  on  the  conductive  solution.  Regions  A,  B,  C, 
D  and  E  are  discussed  in  the  text  (after  Turner,  1973). 
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C,  periodic  solutions  exist  even  though  the  density  stratification  is 
statically  stable.  All  values  in  the  regions  labeled  D  correspond  to 
growth  from  an  unstable  conductive  state.  In  the  upper  right 
quadrant,  or  in  the  diffusive  regime,  values  below  the  Hopf  line 
correspond  to  stable  conductive  solutions  (shaded),  and  points  above 
the  Hopf  line  correspond  to  stable  periodic  ones,  as  denoted  in 
Figure  3.  In  region  E,  where  both  temperature  and  salinity  increase 
with  height,  salt  fingers  are  observed  (Turner,  1973). 

Region  B  in  Figure  4  represents  a  fluid  that  is  gravitationally 
unstable,  since  Ra^,  >  Rag,  but  one  in  which  damped  motion 

nevertheless  does  occur.  This  region  corresponds  to  a 
gravitationally  unstable  case  since  in  a  nondiffusive  fluid  all 
displaced  parcels  will  continue  to  depart  from  their  initial  position 
(Figure  2);  yet  in  this  region  of  Ra^-Rag  space  the  conductive  state 

is  stable.  To  understand  how  this  situation  comes  about,  we  consider 
a  parcel  that  is  perturbed  upwards  in  a  diffusive  fluid.  In  this 
case,  it  is  saltier  and  warmer  than  its  environment,  as  shown  by 
position  1  in  Figure  5.  The  greater,  more  rapidly  diffusing  parcel 
temperature  equilibrates  quickly  with  that  of  the  surroundings  of  the 
parcel  (position  2).  Owing  to  the  salinity  difference  still  present, 
the  parcel  is  heavier  than  its  environment  and  so  acelerates  downward 
toward  position  3B.  The  parcel  will  always  overshoot  its  initial 


Figure  5  Schematic  depiction  of  the  motion  of  parcels  in  the  stable 
conductive  and  periodic  regimes.  A  parcel  is  perturbed 
upwards  from  position  0  to  position  1.  Comparing  the 
parcel  temperature  and  salinity  to  those  of  its  environment 
will  determine  if  the  parcel  returns  to  its  original 
position  via  a  damped  oscillation  (1,2,3A)  or 
(1,2,38,4,5a),  or  maintains  a  periodic  solution 
( 1 , 2 , 3B  4 , 5B ) . 
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position  since  by  Archimedes  Principle  dw/dt  =  g (p  -p  )/p  ,  where  p_ 

e  p  p  p 

is  the  density  of  the  parcel  and  pg  is  the  density  of  the 

environment.  This  acceleration  is  evident  in  (2.9);  since  dw/dt  ' 
<*T’  -PS’,  either  a  negative  temperature  perturbation  or  a  postive 
salinity  perturbation  will  result  in  downward  acceleration.  In 
region  B,  the  salinity  gradient  is  not  large  enough  to  cause  the 
parcel  to  attain  position  3B  in  Figure  5;  it  reaches  only  position 
3A.  This  causes  a  damped  oscillation  to  the  original  position  0. 

In  region  C,  the  fluid  is  stable,  and  Rag  >  Rap.  Again,  a 

parcel  displaced  upwards  is  wanner  and  saltier  than  its  environment 
as  shown  by  position  1  in  Figure  5.  As  before,  the  temperature 
begins  diffusing  rapidly  to  that  of  the  background  state, 
equilibrating  with  its  surroundings  more  rapidly  than  does  the  salt 
field  and  the  parcel  now  moves  to  position  2.  Again,  it  is  heavier 
than  its  environment  and  sinks.  It  travels  below  its  original 
position  to  location  3B,  and  is  now  fresher  and  colder  than  its 
surroundings.  With  the  salinity  gradient  sufficiently  large,  the 
above  process  continues,  and  a  stable  oscillation  is  seen.  Here,  a 
density  distribution  decreasing  upwards  does  not  imply  a  decaying 


solution. 
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2.4  The  Spectral  Model 

To  develop  the  nonlinear,  low-order  spectral  model,  we  need 

functions  for  r* ,  T* ,  and  S*  that  can  describe  the  convective  motion 
with  a  small  number  of  terms  that  satisfy  the  boundary  conditions 
(2.20).  Because  of  the  cyclic,  stress-free  boundary  conditions,  and 
rectangular  domain,  trigonometric  functions  suffice.  With 
wavenumbers  n  =  1  and  n  =  2  in  the  vertical  and  wavenumbers  m,  j , 
m-j,  and  m+j  in  the  horizontal,  the  stream  function  can  be 
represented  in  the  following  equation: 

>*tx*,z*,t*)  =  ^mi(t*)sin(nBc*)sin(z*)  +  ( t* )sin(  jx* )sin( z* )  + 

,  _ ( t* )sin[ (m-j )x* Jsin(2z* )  + 

V_. ,  ( t* )sin{ (m+j )x*)sin( 2z* )  .  (2.50) 

mT  j  ,  2 

The  time-dependent  coefficients  are  amplitudes,  and  these  are  found 
numerically  by  integrating  the  ordinary  differential  system  composing 
the  model.  The  four  harmonics  are  chosen  to  allow  the  possibility 
for  an  exchange  of  stability  between  a  primary  periodic  solution  Pm 

of  wavenumber  m  and  a  primary  periodic  solution  Pj  of  wavenumber  j ; 

this  exchange  occurs  via  the  catalytic  (m-j,2),  (m+j, 2)  terms  at  the 
diffeo-Hopf  points  Dm  and  Dj  (Chang  and  Shirer,  1984;  Wells  and 


Shirer,  1987).  Stability  is  transferred  from  a  solution  dominated  by 
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one  horizontal  wavenumber  to  one  dominated  by  another  wavenunber  and 
can  be  represented  by  a  secondary  branch  emanating  from  a  primary 

one  (Figure  6).  As  occurs  in  rotating  convection,  we  expect  that  a 
two-torus  appears  via  a  secondary  branch  bifurcating  from  one  of  the 
primary  periodic  solutions  (Wells  and  Shirer,  1987). 

From  the  temperature  term  in  the  vertical  velocity  equation 
(2.9)  and  the  vertical  velocity  term  in  the  thermodynamic  equation 

i 

(2.11),  we  see  that  w  ,  or  d^/dx,  and  T  need  the  same  phasing  to 
include  a  heat  flux  term  crucial  to  convection.  Thus,  we  may  write 
the  temperature  field  as 

T* (x* , z* , t* )  =  (t*)cos(mx*)sin(z*)  +  Tj1 (t*)cos( jx*)sin(z*)  + 
T  i  _(t*)cos[(m-j)x*Jsin(2z*)  + 

m-j  t  d. 

T_. H  (t*)cos[ (m+j)x*)sin(2z*)  + 

hit  j  ,2 

Toa(t*)sin(2z*)  +  T04(t*)sin(4z*)  ,  (2.51) 

and  the  salinity  field  as 

S*(x*,z*,t*)  =  Snu(t*)cos(mx*)sin(z*)  +  Sjx (t*)cos( jx*)sin(z*)  + 
Sm_i  O(t*)cos[ (m-j)x*Jsin(2z*)  + 

m—j  t  a 

S_.  ■  <t*)cos[(m+j)x*}sin(2z*)  + 

fflTJ  f  2 

Soa(t*)sin(2z*)  +  S04(t*)sin(4z*)  .  (2.52) 

*  * 

Further,  these  cosine  terms  match  those  from  the  dy  /3x  terms  in  the 
thermodynamic  and  vorticity  equations.  With  representations  this 


Mm  j<l 


JXm'l  ''*»"*»  ♦  Mm,j»l 
Mi*  Mj=l 


Figure  6 


Schematic  bifurcation  diagrams  showing  (a)  periodic  and  toral 
and  (b)  steady  solutions  of  horizontal  wavenumbers  m  and  j. 
The  solid  curves  represent  stable  solutions,  and  the  dashed 
curves  denote  unstable  ones.  The  periodic  branches  Pffi  and  Pj 

are  characterized  by  Hopf  ( )  and  diffeo-Hopf  (D^ )  points 

where  k=m  or  j,  and  a  secondary  toral  solution  S^.  The 

magnitudes  of  the  characteristic  multipliers  Mm  and  Mj 

corresponding  to  the  steady  bifurcation  points  on  the 
conductive  solution  and  the  steady  regular  turning  points  on 
the  steady  solutions  and  Sj  are  shown  in  (b) . 
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general,  J\fw*T*dx*dz*  *  0  is  possible.  Recall  that  background 
profiles  are  horizontal  averages  of  the  variables.  Thus,  it  is  the 
T02,  T04,  S02,  and  SQ4  terms  that  represent  the  nonlinear  distortions 

by  the  rolls  of  the  linear  background  temperature  TQ  and  salinity  S0 

fields.  Also,  the  (m,l),  (j,l),  (m+j,2),  and  (m-j,2)  wavenumber 
pairs  included  in  the  expansions  (2. 50)- (2. 52)  are  only  four  of  the 
infinite  number  of  wavenumber  pairs  represented  by  the  partial 
differential  equations. 

<  ♦  %  • 

Upon  substituting  these  Fourier  series  for  ,  T  ,  and  S  into 
the  dimensionless  equations  of  motion  (2.36)-(2.38) ,  we  can  create 
the  spectral  model  using  the  Galerkin  technique  (Higgins,  1987).  The 
number  of  ordinary  differential  equations  in  this  model  is  determined 
by  the  total  number  of  coefficients  in  the  expansions — sixteen  in 
this  case: 


Smi=  MunSo*  -  4(m-J)*USm+j,2  '  4<m+j)Vj,aSjl  +  4<m+J >*jiSm-j ,a  " 
4(m-j>Vi,aSjl+  -  ismi(m2a2+1)  ■  (2*53) 


Sjx=  j*jxS02  +  4(m+j>VSm-j,a  +  4lm-j)Sm+j,aV  +  4(m+j  )Snu  Vj  ,a  + 


4<m'j)^m+j,aSmi+  -  iSji(J2a2+1)  * 


(2.54) 
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+  4]  , 


"  4lm  +  4(m-J)>'j1SII11  +  2(“+J^m+j>2S04  + 

~  ISm+J(2l<ma+Ja)2  +  4]  , 

S02=  ~  2  “  2  & j ^  ^  SQ2  , 


So*= 


m+j,2bm«fj,2  "  16aS 


04 


4(m-J^m+Jf2TJi+  mR^rmi  -  ^Tmi(m2a2+1)  , 


+  |(m+j)Vji 


^jx  J>jiToa  +  4<m+J)rmiTm-j,a  +  4l®-j,Tm+j,a*mi  +  4<m+j>Tmi 

jRVj,  -  ^Tjl { j*a*+l )  , 

Vj,s-  -  'll,"*jl,raiTji  ■  4l»*-JlV-jlTnli  ♦  2,“‘J  1  Vj ,2T04  * 
<m-JISsTVj,2  -  4  4]  , 


T  =  _  1 
®+J»a 


4(®-J^m1Tj1  +  4(n>-J^j1Tmi  +  2<m+j>T('m+j>2T04  + 
(m+jjRaj,^  a  -  aTm+j,2f(na+ja>2  +  4]  , 


(2.55) 


(2.56) 

(2.57) 

(2.58) 

r 

m-j.a 

(2.59) 

®-j.2  + 

(2.60) 


(2.61) 


(2.62) 


I 
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^02"  m  2  ^miTmi  “  2  ~  a  ^02  ’ 


(2.63) 


T04=  ~  <m+j>Vj,aVj,2  “  16a  To«  •  (2'64) 


;  =  ^!ffl^(m3a2_m2Ja2_2lIlJ2a2+3m+3J)  +  rAllB±J.  1 2  j  -m3a2-  m*a2j+ 

mi  4(m2a2+l )  4(m2a2+l ) 


2mjaa2-3m+3j)  +  -  ^<m2a2+l)  , 

In  a  +1 ) 


(2.65) 


v  y 

LnuJ" 


Vi.=  -E^rJBr*L,A(-j3a2+inj2a2+2jm2a2-3m-3j)  +  ■B-1-Bl^J-t2(-j3a2-  j2a2m+ 


J1  4( j2aa+l ) 


4( j2a2+l ) 


P*  : 


2 jm2a2+3m-3 j )  +  — - (T.-S-  ) - ^(j^+l)  ,  (2.66) 

( j2a2+l )  Ji  J1  a 


•  _  t  PNKV.Al,) 

m  j’2  4 [ a2 ( m- j ) 2  +  4]  a2(m-j)2  +  4 


-  I  Vj,2ta2(m-j)2  +  41  * 


(2.67) 


.  _  j-aj  +j  i  ,  Pi->jnyj.,-s^j.,> 

ni+j.a  4(aa(m+j)a+  4]  aa(nH-j)a  +  4 


a  ^m+j,2[aa(m+J,a  +  41  ' 


(2.68) 


The  modification  of  the  background  temperature  and  salinity 
fields  arises  from  the  interaction  of  modes  via  the  nonlinear,  or 


Jacobian  terms,  in  the  equations  of  motion.  For  example,  the 
nonlinear  product  TJniV'mi  modifies  the  Toa  term  in  (2.63),  since 
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adding  the  vertical  wavenumbers  and  subtracting  the  horizon k 
wavenumbers  gives  the  02  mode.  Likewise,  for  TQ4,  the  product 

Tm+j ,2^m+j ,2  niodifies  the  TQ4  term  in  (2.64).  A  higher  order 

spectral  model  would  be  necessary  to  produce  changes  in  terms  such  as 
TQ1  and  TQ3 .  Thus,  a  nonlinear  model  is  needed  to  represent  the 

changes  in  the  mean  background  temperature  and  salinity  profiles. 

The  above  differential  system  (2. 53)- (2. 68)  contains  two  primary 

periodic  solutions  of  the  form  in  the  Veronis  (1965)  five-coefficient 

model  V  , 
m 


Pjd  -  , 0 , 0 , 0 ,  SQ2 , 0 >Tmt , 0 , 0 , 0 ,  T02  ^ )  >  (2.69) 


and  that  in  V • , 


Pj  =  (0,SJl,0,0,So2,0,0,TJl,0,0,Toa,0,0,^Jj,0,0),  (2.70) 

as  well  as  two  other  periodic  solutions 

P»,j:(0-0'Vj,a*0-0-Sod0-0-Vj,2-0'0*T04-0-0-Vj,a-0>  ’  ,2-71) 

Bind 

W0’0’0’SWO,S’0-So.’0'0'0’Vj.S-0-T..'0'0-0-Vj.a»  •  <2-72) 

The  five-coefficient  steady  solutions  from  the  Veronis  model  su'e  also 
included  in  this  truncation. 


In  the  Veronis  model,  the  primary  branches  (2.69)  and  (2.70) 


bifurcate  from  the  conductive  solution  at  the  Hopf  points  that  were 
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calculated  using  (2.47)  for  the  partial  differential  equations.  The 

five-coefficient  periodic  solutions  P„  and  P-  will  be  essential  for 

m  j 

deciding  whether  toral  or  other  solutions  are  possible.  We  will  find 
these  solutions  using  an  eigenvalue  analysis,  with  which  we  may  also 
determine  the  stability  of  the  periodic  solutions  Pm  and  Pj .  With 

this  stability  analysis,  we  can  show  when  either  five-coefficient 
branch  loses  stability  to  perturbations  of  different  wavenumbers ,  as 
represented  by  the  other  eleven  coefficients.  From  this  eigenvalue 
analysis,  we  can  also  diagnose  what  type  of  solution  gains  stability 
when  the  five-coefficient  branch  loses  stability.  Of  particular 
interest  would  be  the  toral  solution  that  gains  stability  at  a 
dif feo-Hopf  bifurcation.  These  properties  are  discussed  in  Chapter 


3. 
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Chapter  3 

FINDING  BIFURCATION  POINTS  ON  PERIODIC  SOLUTIONS 

To  discover  when  a  new  solution  bifurcates  from  a  temporally 
periodic  one,  we  need  to  determine  the  stability  of  the  temporally 
periodic  branch  (Feldstein,  1987).  It  is  sufficient  to  show  that  a 
new  solution  is  created  when  this  primary  branch  loses  stability  at  a 
certain  value  of  the  forcing  rate.  Determination  of  both  the 
stability  of  the  periodic  solution  and  the  character  of  the  new 
branching  solution  is  performed  via  an  eigenvalue  analysis  of  an 
appropriate  system  of  perturbation  equations.  By  examining  these 
eigenvalues,  or  characteristic  multipliers  ^n,  as  the  forcing  rate  is 

varied,  we  are  able  to  find  diffeo-Hopf  points  that  signal 
bifurcation  to  a  torus,  as  well  els  other  bifurcation  points,  on  the 
periodic  solution.  The  behavior  of  the  characteristic  multipliers  at 
the  bifurcation  points  reveals  the  type  of  new  solution  emanating 
from  the  periodic  one. 

3.1  Stability  of  a  Periodic  Solution 

To  determine  the  stability  of  a  periodic  solution,  we 
superimpose  a  perturbation  on  that  solution  and  determine  whether  the 


perturbation  grows  or  decays.  If  this  perturbation  grows,  then  the 
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periodic  solution  is  unstable,  and  if  it  decays,  then  the  solution  is 
stable . 

We  can  write  the  sixteen-coefficient  model  (2. 53 >-(2.68)  in 
the  form 

yi  =  fi(y1  ,y£ , . . .  ,y16  (RagjRap^.r.P,!!!,  j )  ,  i=l,...,16  ,  (3.1) 

which  has  a  periodic  solution  y^(t)  =  P^(t)  of  period  T  (Feldstein, 
1987).  We  then  consider  perturbations  x^(t)  of  P^(t)  and  write 

y^(t)  =  Pi ( t )  +  xi(t)  ,  i=l . 16  .  (3.2) 

By  knowing  whether  | x^ ( t ) |  increases  or  decreases ,  we  can  determine 
the  stability  of  P^lt). 

To  obtain  a  set  of  equations  that  describes  how  these 
perturbations  change  with  time,  we  substitute  (3.2)  into  (3.1).  For 

P  P 

example,  we  let  Smi(t)  =  Smi  +  xj5  tmi(t)  =  Vmi  +  xJ3,  etc.,  where 

the  superscript  p  denotes  a  periodic  solution,  and  substitute  these 
into  (2.53).  To  linearize  about  P^(t),  we  assume  that  |x^(t)| 

remains  small  and  so  drop  products  of  the  perturbations.  After 
equating  the  periodic  and  perturbation  terms,  we  get 


P  P 


Smi(t)  =  -  +  4(m+j,Sm-j,2*jx 


.  -1( 


4(m-j)Vj,aSji  +  “VSq*  +  m  "vC  "  <T/a>S^<n,aa8  +  D  ,  <3. 3) 
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and 


x1(t)  = 


-  ?<m-J>[xi4S^J>2  +  x4^]  -  i(m+j)[xi3S^  +  v£-j,a]  + 


4(m+j)[xi4Sm-j,2  +  Vji]  -  J,«-j)[vLjf2  +  Xi6SJij  +  mRasxi: 


P  P 

+  «“(X13S02  +  x8vmi)  -  (r/a)xi (m2a2+l)  .  (3.4) 

We  note  that  (3.3)  is  the  periodic  version  of  (2.53),  and  that  the 
last  three  terms  of  (3.3)  are  those  of  the  Veronis  model  Vm.  Also, 

(3.4)  is  linear  in  the  perturbation  variables  x^.  The  other  fifteen 

pairs  of  equations  are  formed  in  a  similar  way. 

The  resulting  thirty-two-coefficient  system  S  consists  of  the 
original  sixteen-coefficient  model  (2.53)  -  (2.68)  and  a  set  of 
sixteen  perturbation  equations.  The  system  of  perturbation  equations 
admits  of  the  trivial  solution,  while  the  first  sixteen  give  the 
periodic  solution.  To  find  the  behavior  of  the  perturbation 
components  x^,  we  need  not  know  the  periodic  solution  analytically. 

In  this  case,  the  first  sixteen  equations  give  the  periodic  solution 
numerically,  and  these  values  feed  into  the  sixteen  perturbation 
equations  for  the  stability  calculation.  By  calculating  the 
eigenvalues  pn,  n=l,...,16,  of  a  matrix  derived  from  solutions  to 

this  system  evaluated  at  one  period,  we  can  determine  if  the 
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perturbation  of  the  periodic  solution  in  the  sixteen-coefficient 
model  is  growing  or  decaying  in  any  of  the  sixteen  dimensions 
(Feldstein,  1987).  We  outline  how  we  calculate  the  values  of  these 
characteristic  multipliers  in  the  next  section. 

3.2  Eigenvalue  Calculation  and  Analysis 

Given  a  five-coefficient  periodic  solution  Pm  or  Pj ,  we  would 

like  to  know  those  parameter  values  for  which  this  solution  is  stable 
within  the  sixteen-coefficient  model,  those  for  which  it  is  unstable, 
and  the  type  of  solution  that  replaces  the  periodic  one  when  it 
becomes  unstable.  From  the  eigenvalue  analysis  of  the  integrated 
perturbation  system  discussed  in  Section  3.1,  all  of  this  information 
can  be  determined. 

The  initial  conditions  for  each  integration  of  the  thirty- two- 
equation  system  S  e  e  specified  in  the  following  way.  The  first 
sixteen  coefficients  of  S  are  chosen  to  be  on  the  primary  branch  Pm 

or  P  which  is  being  studied.  Because  these  branches  are  solutions 
J 

to  the  full  system,  the  appropriate  five  Veronis  coefficients  are 
nonzero,  the  other  eleven  coefficients  are  zero,  and  these  eleven 
will  not  grow  numerically.  The  equations  composing  the  second  half 
of  S  are  set  to  the  value  one  in  the  perturbation  component  x^  with 
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all  others  set  at  zero;  the  result  is  integrated  for  one  period  T. 
This  procedure  is  repeated  for  each  value  of  k=l , 16  to  produce 
final  values  [x^t)]^,  i=l,...,16,  after  one  period.  These  values 

are  written  as  the  entries  of  a  matrix  B.  For  example,  to  find  the 
first  column  k=l,  the  initial  conditions  for  the  perturbation 

T  T 

equations  are  (x. ,x0 , . . . ,x„ . )  =  (1,0, ...,0)  .  These  integrations 

1  2  X  0 

are  performed  sixteen  times,  one  for  each  set  of  initial  conditions 
involving  only  one  nonzero  perturbation  component .  The  resulting 
sixteen  columns  yield  the  matrix  B  of  perturbations  evaluated  at  one 
period  T: 


:  1 

2 

16 

b 

b  „ 

.  .  .  b 

li 

12 

1  f  i  0 

b 

b  „ 

.  .  .  b 

21 

22 

2,10 

• 

# 

(3.5) 


*-  bie,i  bie,«  *  '  *  bi 6,i6 

The  eigenvalues  Pn  of  matrix  B  give  the  desired  stability 

information.  Recall  that  if  these  eigenvalues  all  have  magnitude 
less  than  one,  then  the  periodic  solution  is  stable  (Feldstein, 
1987).  As  we  note  below,  this  matrix  can  be  separated  into  a  5x5,  a 
9x9,  and  two  lxl  matrices  without  changing  the  values  of  the 
characteristic  multipliers.  Oie  eigenvalue  of  B  will  equal  one, 
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indicating  integration  for  one  period.  If  the  model  is  integrated 
for  one  time step  fit  instead  of  one  period,  then  the  eigenvalues  of 
matrix  B(fit)  give  instantaneous  divergence  rates  (Nese,  1989). 

After  determining  the  Hopf  points  and  Hj  analytically  from 

(2.47),  we  can  integrate  the  five-coefficient  Veronis  model  Vffl  or  Vj 
for  values  of  Ra^,  that  are  slightly  greater  than  those  for  or  ttj 
to  f ind  a  periodic  branch  Pffl  or  P j  ( cf .  Figure  6 ) .  Because  the  f i ve 
nonzero  components  of  the  Veronis  model  V^,  k  =  m  or  j,  also  produce 

a  solution  to  the  sixteen-coefficient  model  ( 2.53 )-(2.68)  with  the 
other  eleven  components  set  to  zero,  it  is  much  easier  and  faster  to 
integrate  the  five-coefficient  model  than  the  original 

sixteen-coefficient  model.  Whenever  we  are  able  to  find  Pm  or  Pj 

numerically,  then  it  follows  that  these  solutions  are  stable  relative 
to  perturbations  involving  the  five  components  of  the  Veronis  model. 

After  integrating  the  sixteen-coefficient  system,  we  may 
separate  it  into  four  subsystems  linearized  about  Pffi  or  Pj .  That  is, 

we  may  rewrite  the  linearized  system  as 


c 

c 

0 

c 

d 

0 

D 

0 

d 

• 

V  *  ^ — 

x  dt 

♦ 

-*  o 

GO 

= 

0 

-16-r/a 

0 

S04 

» 

» 

T 

04 

0 

-16/a 

T 

04 

(3.6) 


Here,  c  denotes  perturbations  of  the  components  of  V^,  d  denotes 

i  t 

perturbations  in  the  9x9  system,  and  SQ4  and  TQ4  represent 
perturbations  of  the  S04  and  T04  terms,  respectively.  From  (3.6)  it 

)  I 

is  easy  to  see  that  SQ4  and  TQ4  decay  to  zero  according  to 

So*  =  S^4(0)e(_16rt/a)  ,  (3.7) 

and 

T04  =  T„4(0)e(-16t/a)  .  (3.8) 

Spectrally,  the  periodic  solutions  will  always  be  stable  to 
perturbations  in  these  directions;  the  corresponding  eigenvalues  will 
always  have  magnitude  less  than  one.  Now  the  stability  analysis  has 
been  reduced  to  evaluating  the  properties  of  the  matrices  C  and  D. 
Here,  C  is  a  5x5  time-dependent  matrix  containing  the  terms  from  Vm 

or  V  • ,  and  D  is  a  9x9  matrix.  As  a  check,  C  and  D  can  also  be  found 

J 

by  separate  integrations,  and  the  eigenvalues  of  these  matrices 
should  be  the  same  as  those  found  when  integrating  the  entire  system 
S.  Since  the  single  wavenumber  models  Vm  or  Vj  cannot  contain  a 

toral  solution,  the  stability  information  for  secondary  branching  to 
a  torus  is  found  in  matrix  D. 

Since  the  eigenvalues  measure  expansion  or  contraction  relative 
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to  the  orbit  Pm  or  Pj ,  one  eigenvalue  ,  which  measures  variations 

in  the  direction  of  the  orbit,  should  always  equal  one;  this 
eigenvalue  is  found  in  the  5x5  submatrix.  The  existence  of  p  =1 

provides  a  very  useful  check  on  the  accuracy  of  the  calculations.  If 
is  slightly  greater  than  or  less  than  one,  then  it  may  indicate 

that  the  timestep  used  in  the  integration  is  too  large  or  the 
number  of  iterations  is  too  small,  causing  the  integration  to  stop 
just  before  or  just  after  one  period. 

If  the  periodic  solution  is  determined  numerically,  then  the 
period  can  only  be  known  approximately.  Therefore,  a  careful 
analysis  of  the  magnitude  of  the  period,  and  the  appropriate  number  N 
of  intervals  needed  to  divide  it,  is  necessary  to  determine  the 
optimal  timestep  At.  Here,  the  number  of  required  timesteps  is 
usually  between  13,000  and  20,000,  with  a  timestep  At  on  the  order  of 

10-4.  This  is  sufficient  to  yield  =  liO.Ol  and  to  approximate 

|AVkll  =  [^ki(t=0)  "  (t=0)  10  order  10”4*  The 

sensitivity  of  the  five  largest  eigenvalues  is  illustrated  in  Table 

1. 

Table  1  contains  eigenvalues  calculated  using  two  different 
timesteps  At  and  lengths  of  integration  N.  The  best  approximation  is 
the  one  in  the  third  colum,  since  is  the  smallest.  Here,  the 
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Table  1  Estimates  of  the  five  largest  eigenvalues  at  for  Rag  =  5.0, 
Ra^,  =  12.8,  ^  =  12.25,  a  =  0.486,  m  =  2,  and  j  =  1  for 

different  timesteps  At  and  iterates  N  in  the  integration. 
The  first  three  eigenvalues  are  from  integrating  the  5x5 
subsystem  involving  matrix  C  in  (3.6)  for  the  Veronis  Pm 

solution,  and  the  next  three  are  from  the  9x9  subsystem 
involving  the  matrix  D.  Also,  T„  is  the  approximate  period 

a 

where  T  =  NAt  (cf.  Figure  8). 

8 


At 

3xl0-4 

3xlQ“4 

2.2xl0~4 

2.2xl0“4 

N 

15366 

15360 

20953 

20954 

Ta 

4.6098 

4.6080 

4.6096 

4.6099 

|AVkil 

2.40xl0-3 

3. 15xl0-2 

2.15xl0-4 

3.89xl0“3 

A> 

1.01960 

0.99368 

1.0179 

1.0205 

1 

-0. 025371 

P 

0.96993 

0.99368 

0.97146 

0.96911 

2 

+0.025371 

3 

0.09332 

0.09323 

0.09331 

0.09333 

A* 

4 

2.01620 

2.01100 

2.01580 

2.01640 

-0.09066 

-0.089793 

-0.090600 

-0.09070 

±0. 324141 

±0.324731 

±0.324191 

±0.324111 
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eigenvalue  that  indicates  one  period  of  integration  is  .  After 

integrating  one  more  times tep  (colunn  four),  we  see  that  j I  has 

increased  15.5  times,  indicating  a  loss  of  two  decimal  places  of 
accuracy  in  the  estimate  of  (T).  Also,  the  most  important 

eigenvalues,  which  are  those  close  to  one  in  magnitude,  change  in  the 
second  decimal  place  as  the  model  is  integrated  an  extra  times tep. 

We  note  that  v  must  be  known  to  four  decimal  places  to  get  the 
eigenvalues  accurate  to  two  decimal  places. 

As  the  timestep  is  increased  to  3.0xl0~ 4  ( columns  one  and  two), 
the  behavior  of  the  first  two  eigenvalues  is  deceiving.  As  the 
number  of  integrations  decreases,  the  first  two  eigenvalues  become 
complex  and  have  a  magnitude  near  one,  which  is  indicative  of  a 
diffeo-Hopf  bifurcation.  Knowing  that  we  must  have  one  eigenvalue 
equal  to  one,  however,  we  conclude  that  this  cannot  be  a  complex 
conjugate  pair.  Also,  the  imaginary  part  is  quite  small.  With  a 
better  approximation  to  the  period,  we  see  that  one  of  these 
eigenvalues  is  signaling  the  occurrence  of  either  a  harmonic 
bifurcation  or  a  periodic  regular  turning  point,  since  we  are  far 
from  the  Hopf  bifurcation  point,  and  the  other  indicates  that  we  have 
integrated  for  one  period.  This  conclusion  would  change  if  we  were 
close  to  H^.,  because  in  this  case  there  is  also  an  eigenvalue  having 
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magnitude  near  one  and  decreasing  as  we  move  away  from  The 

errors  in  the  values  of  that  are  illustrated  in  Table  1  become 
greatly  magnified  for  large  values  of  lAY^I,  i.e.  before  an  accurate 

value  of  the  period  has  been  determined. 

The  other  fifteen  eigenvalues  determine  the  stability  of  the 
primary  branch,  and  depend  sensitively  on  how  accurately  the  value  of 
the  period  is  known.  These  characteristic  multipliers  also  indicate 
which  type  of  new  solution  has  bifurcated  from  the  primary  branch. 
Three  types  of  new  solutions  are  possible,  as  we  review  below. 

3.2.1  Subhannonic  Bifurcation 

If  a  real  eigenvalue  pt=-l  occurs,  then  a  subharmonic ,  or  period 
doubling,  bifurcation  is  signaled  (Ruelle,  1989).  The  eigenvalue 
=— 1  implies  that  an  eigenvector  of  the  linearized  system  has  been 

rotated  180°  after  one  period  T.  Thus,  the  orbit  has  twice  as  far  to 
go  to  complete  one  cycle.  Because  it  requires  twice  the  time  to 
return  to  its  original  position,  the  bifurcation  is  called  a  period 
doubling  bifurcation  and  the  period  of  the  new  solution  is 
approximately  2T.  In  this  case,  there  is  only  one  new  periodic 
solution  that  gains  stability  from  the  old  one  (Sparrow,  1982; 


Feldstein,  1987). 
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3.2.2  Harmonic  Bifurcation 

A  real  eigenvalue  of  indicates  a  harmonic  bifurcation 

(Ruelle,  1989).  This  eigenvalue  represents  an  eigenvector  that 

rotates  360°  in  one  cycle  and  therefore  does  not  change;  this  implies 
that  the  new  periodic  solution  will  have  the  same  period  as  the  old 
one.  In  fact,  t^o  new  periodic  solutions  bifurcate  from  the  old 
solution  (Sparrow,  1982;  Feldstein,  1987). 

3.2.3  Toral  Bifurcation 

A  toral  or  diffeo-Hopf  bifurcation  occurs  when  a  complex 
conjugate  pair  of  eigenvalues  has  magnitude  one.  Here,  Mg  = 

exp(i2TT6),  0  <  9  5  1,  and  9  can  be  rational  or  irrational.  If  9  is 
rational  and  9  *  1,1/2, 1/3,  or  1/4  (Iooss  and  Joseph,  1980; 
Ouckenheimer  and  Holmes,  1983),  then  the  resulting  solution  is  a 
periodic  one  residing  on  the  surface  of  a  two- torus.  If  9  is 
irrational,  then  the  orbit  covers  the  surface  of  a  two-torus  and  is 
called  quasi -periodic. 

Quasi -periodic  solutions  or  tori  are  expected  to  appear  in  many 
routes  to  chaos  (Nese,  1987).  As  shown  by  Wells  and  Shirer  (1987), 
tori  can  connect  two  periodic  solutions,  such  as  Pm  and  P^,  of 


50 


different  horizontal  wavenumber.  The  corresponding  bifurcation 
diagram  takes  on  a  goal  form  as  in  steady  convection  (Hirschberg, 
1987).  Therefore,  exciting  two  wavenumbers  one  after  the  other  by 
increasing  the  forcing  rate  in  a  fluid  may  lead  to  toral  solutions 
(Figure  6a) . 

The  above  bifurcations  are  straightforward  to  track  by  following 
the  changes  in  magnitudes  of  the  eigenvalues.  Even  unstable  periodic 
solutions  in  the  sixteen-coefficient  model  can  be  traced  if  the 
instability  originates  from  an  eigenvalue  of  the  9x9  submatrix  D  in 
(3.6).  In  the  next  chapter,  we  will  use  these  ideas  to  investigate 
whether  toral  bifurcations  can  occur  in  thermohaline  convection. 
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Chapter  4 

PHYSICAL  INTERPRETATION  OF  MODEL  SOLUTIONS 

Vv'e  begin  integrating  the  model  for  values  of  Ra^,  >  as 

discussed  in  Chapter  2.  The  results  of  these  integrations  give 
stable  steady  or  periodic  solutions.  By  examining  the  stability  of 
these  solutions,  we  can  determine  when  toral  solutions  bifurcate  from 
the  periodic  ones  and  whether  these  new  toral  solutions  are  stable  or 
unstable.  In  this  chapter,  specific  solutions  for  realistic 
parameter  values  are  examined  to  determine  the  changes  in  density, 
temperature  and  salinity  characterizing  these  solutions.  In 
particular,  we  focus  on  the  changes  in  the  background  profiles  of 
density. 

4 . 1  Bifurcation  Diagrams  Representing  Model  Solutions 

To  construct  bifurcation  diagrams  for  particular  values  of  Rag 

and  a  that  show  the  steady  and  periodic  solutions  of  the  model,  we 
find  the  periodic  branches  numerically,  while  we  determine  the  steady 
solutions  analytically.  To  find  the  forcing  rates  yielding  toral 
solutions,  we  begin  integrating  the  model  using  values  of  Rag  and  Ra^- 

that  are  just  past  those  tor  the  bifurcation  points  1^  and  H^;  these 

integrations  give  the  Pm  and  Pj  branches,  as  shown  in  Figure  6a. 


52 


Using  the  eigenvalue  analysis  discussed  in  Chapter  3,  we  look  for 
diffeo-Hopf  bifurcation  points  on  both  the  Pm  and  Pj  branches;  as  we 

noted  in  Section  3.2.3,  these  bifurcation  points  signal  the  existence 
of  toral  branches. 

Diffeo-Hopf  points  emanate  from  the  double  Hopf  point  on  the 
conductive  solution  as  Pffi  and  Pj  are  separated.  These  points  are 

denoted  D,  on  Figures  7  and  8.  Since  the  steady  convective  solutions 

Sm  and  Sj  generally  have  large  basins  of  attraction,  we  would  like  to 

move  the  periodic  branches  and  diffeo-Hopf  points  as  close  as 
possible  to  the  steady  regular  turning  points  Fm  and  Fj ,  where  the 

basins  of  attraction  for  and  S.  are  the  smallest  and  the  basin  for 

m  J 

the  toral  solution  is  the  largest,  provided  the  toral  solution  ip 
stable . 

Reducing  the  basins  of  attraction  for  the  steady  solutions 
and  Sj  can  be  achieved  by  making  the  value  of  less  than  that  of 
both  Fm  and  Fj,  which  are  the  minimum  values  of  Ra<p  for  the  steady 
solutions  to  the  models  and  Vj,  respectively  (Figures  7  and  8). 

The  equation  describing  the  steady  states  can  be  calculated  by 
setting  the  temporal  derivatives  of  the  five-component  model,  k  = 
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Figure  7  Bifurcation  diagram  for  Rag=1.6836,  a=0.486,  m=2,  j=l , 

P=7  and  r=0.01.  Hopf  and  diffeo-Hopf  bifurcation  points 
are  denoted  by  and  Dj ,  respectively,  where  k=m  or  j. 

Here,  and  Gjg  represent  periodic  regular  turning 

points  on  the  periodic  branches  PB  and  P^,  and  Fm  and  Fj 

are  steady  regular  turning  points  on  the  steady  branches 
Sm  and  Sj.  There  is  no  stable  secondary  branch  S^. 

Solid  curves  represent  stable  solutions,  and  dashed 
curves  denote  unstable  ones. 
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Figure  8  Bifurcation  diagram  for  Rag=5.0,  a=0.486,  m=2,  j=l,  P=7 

and  rrO . 01 .  Hopf  and  diffeo-Hopf  bifurcation  points  are 
denoted  by  and  Dj,  respectively,  where  k=m  or  j. 

Here,  G^g  and  represent  periodic  regular  turning 

points  cm  the  periodic  branches  P_  and  P  • ,  and  F  and  F 

m  J  ID  J 

are  steady  regular  turning  points  on  the  steady  branches 

and  S  ^ .  There  is  no  stable  secondary  branch  S  - . 
v  mj 

New  bifurcation  points  are  the  harmonic  ones  G^  and  Gjg 
on  Pm  and  Pj.  Solid  curves  represent  stable  solutions, 
and  dashed  curves  denote  unstable  ones. 
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m  or  j ,  to  zero  and  then  by  manipulating  the  remaining  terms  to 
produce  a  quadratic  in  <V 


fjlikVt*  +  T^Js^akr2  +  8f  3ak  -  8a3RaTk3  +  8a3k3rRag]  + 

64*4r2/ak  -  64aRapkr2f  +  64TaRagk£  =  0  ,  (4.1) 


where  $sk*aa+l.  This  equation  gives  the  parabolic  curves  that 
represent  the  m  and  j  steady  solutions ,  Sffl  and  Sj  ( Figures  7  and  8 ) . 

To  obtain  the  value  of  RaT  at  F^,  we  make  the  discriminant  of  (4.1) 
vanish  and  solve  for  Ra^lF^): 


RaT(Fk) 


(rRag)172  + 


( 1-r2 )  Q^a2*!  )3 


k^a2 


1/2 


(4.2) 


Requiring  that  IRa^fH^)!  <  IR&pfF^H  for  both  k  =  m  and  k  =  j 
produces  a  quadratic  inequality  in  Rag  having  a  double  root.  It  is 


interesting  to  note  that  using  the  approximate  formula  (2.48)  for 


instead  of  the  exact  one  (2.47)  gives  poor  results.  Therefore,  use 
of  the  full  Hopf  equation  (2.47)  is  necessary  to  produce  the  correct 
values.  For  m=2,  J=l,  a=  a^  =  0.49343,  and  Rag  =  1.6836,  we 


find  that  H^,  Fm,  and  Fj  all  coincide;  however,  RaT(Hmj)  <  FB  or  Fj 

is  not  possible  for  any  values  of  the  parameters.  Thus,  we  conclude 
that  Rag  =  1.6836  is  the  best  value  of  the  forcing  rate  that  yields 
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values  of  and  Hj  that  are  as  close  to  Fm  and  Fj  as  possible  when 
the  value  of  a  is  varied  from  a^j.  When  the  magnitude  of  the  aspect 
ratio  is  chosen  to  be  slightly  smaller  than  that  of  (as  in 
Figures  7  and  8 ,  where  a  =  0 . 486 ) ,  separates  into  1^  and  Hj ,  and 
these  points  move  to  larger  values  of  R&p.  Also,  Fm  and  Fj  separate 
so  that  Ra^(Fffl)  <  RaT(Fj)  for  a  <  a^.  The  result  is  that  the  steady 
regular  turning  points  always  occur  at  smaller  values  of  Ra^  than  the 
Hopf  points  when  a  *  amj‘ 

The  dashed  lower  branches  of  the  steady  curves  in  Figures  7  and 
8  are  unstable.  However,  from  an  eigenvalue  analysis  of  the  steady 
solutions,  we  find  that  both  upper  branches  are  stable  in  the 
sixteen-coefficient  model.  For  the  branch  S  •  to  be  stable,  even  when 

4/ 

it  bifurcates  from  the  conductive  solution  at  larger  values  of  Ra^. 

than  does  Sm,  there  must  be  a  bifurcation  point  on  the  unstable 

branch  of  the  j  steady  solution  (e.g.  Laufersweiler,  1987; 

Hirschberg,  1987).  Here,  an  eigenvalue  is  equal  to  one. 

This  deduction  can  be  seen  by  tracing  the  eigenvalues  on  the 
steady  solutions  (Figure  6b).  We  note  that  the  magnitudes  of  all  the 
characteristic  multipliers  will  be  less  than  one  on  a  stable  steady 


solution,  but  at  least  one  multiplier  has  magnitude  greater  than  one 
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on  an  unstable  solution;  moreover,  these  magnitudes  vary 
continuously.  To  trace  these  eigenvalues,  we  consider  a  case  without 
Hopf  bifurcations  and  begin  on  the  stable  conductive  solution  before 
the  steady  bifurcation  points  (Figure  6b).  Here,  all  eigenvalues 
have  magnitude  less  than  one,  since  both  the  m  and  j  conductive 
solutions  are  stable.  Then,  we  follow  a  particular  steady  solution, 
m  or  j ,  noting  that  an  eigenvalue  will  have  magnitude  one  at  a 
bifurcation  point  and  the  steady  regular  turning  points  (cf. 

Figures  6b,  7  and  8). 

Since  both  upper  branches  are  stable  in  the  sixteen-coefficient 
mode.  ,  all  sixteen  eigenvalues  have  magnitude  less  than  one  on  these 
solutions.  This  is  found  by  treating  the  steady  solutions  as 
periodic  ones,  and  performing  an  eigenvalue  analysis.  It  is  not 
possble  for  the  S 1  branch  to  be  stable  unless  the  eigenvalue  Mm, 

whic1  became  greater  than  one  at  the  m  steady  bifurcation  point, 

becc  es  less  than  one  before  F:.  Therefore,  there  must  be  a 

J 

bifurcation  point  somewhere  on  the  unstable  branch  of  Sj  at  which 
Pm=i.  This  secondary  steady  bifurcation  may  produce  an  unstable 

secondary  branch  connecting  the  m  and  j  steady  solutions;  as  noted  by 
Chang  and  Shirer  (1984),  such  connecting  branches  are  comnonly  found 


in  convective  flows. 
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As  we  can  see  from  Figure  7 ,  there  are  several  striking 
characteristics  of  the  periodic  solutions.  First,  as  signalled  by 
the  occurrence  of  an  eigenvalue  equal  to  one  in  the  Vm  and  Vj  models, 

the  Pffl  and  Pj  branches  end  suddenly  via  periodic  regular  turning 
points  k  =  m  or  j  (Figure  9a).  This  bifurcation  is  not  a  stable 

harmonic  one  since  no  periodic  solution  replacing  Pffi  can  be  found  if 
the  model  is  integrated  for  values  of  Ra^  >  G^0 ;  the  preferred 
solution  becomes  the  steady  one,  S^,  depending  on  the  periodic  branch 

near  which  the  model  was  initialized. 

The  second  important  characteristic  of  Figure  7  is  the  absence 
of  a  connecting  secondary  branch  indicating  the  existence  of  a  stable 
toral  solution.  A  diffeo-Hopf  point  Dj  is  found  only  on  the  Pj 

branch;  at  this  point  the  periodic  solution  gains  stability.  This 
toral  bifurcation  is  important  since  in  its  absence,  there  would  be 
no  stable  P  solution.  Near  D^,  any  energy  put  into  the  (m,l), 

( m+ j ,  2 ) ,  ( m- j , 2 ) ,  and  (0,4)  wavenunbers  quickly  decays .  Since  no  Dm 

bifurcation  point  exists  on  Pm,  there  can  be  no  connecting  branch  to 

provide  an  exchange  of  stability  between  periodic  solutions  of 
wavenumber  m  and  those  of  wavenunber  j. 

As  the  forcing  rate  owing  to  salinity  is  varied  from  Rag  = 
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(b)  A 


Figure  9  A  schematic  depiction  of  the  periodic  branch  Pm  in  (a) 

the  five-coefficient  model  and  in  lb)  the  sixteen- 
coefficient  modei.  At  the  periodic  regular  turning 
point  G^,  the  stabie  periodic  solution  i solid)  turns 

towards  smaller  values  of  Rar-  and  oecomes  unstable 

(dashed).  At  the  harmonic  bifurcation  point  Gmy ,  the 

periodic  solution  becomes  unstable  owing  to  excitement 
of  some  of  the  other  11  modes. 
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1.6836  to  Rag  =  5.0  (Figures  7  and  8),  the  Hopf  bifurcation  points 
move  to  larger  values  of  Rap,  and  all  solutions  extend  to  larger 
values  of  ( I .  Characteristics  similar  to  those  seen  in  Figure  7 
are  found  at  this  larger  value  of  Rag,  with  the  addition  of  harmonic 
bifurcation  points  C^g  and  Gjg  on  Pm  and  (Figures  8  and  9b).  From 
an  eigenvalue  analysis,  we  find  that  near  Rap  =  12.65  and  Rap  = 

12.75,  an  eigenvalue  from  each  of  the  9x9  m  and  j  submatrices  is  real 
and  has  magnitude  one.  Numerical  evidence  supporting  the  occurrence 
of  a  harmonic  bifurcation  near  the  end  of  Pffi  for  Rag  =  5.0  is  shown 

in  Table  1  in  Chapter  3.  Near  the  end  of  Pffi  at  Rap  =  12.8,  there  is 
an  eigenvalue  approaching  one,  namely  n .  The  other  multiplier 

having  a  value  near  one  indicates  integration  for  one  period.  We 
recall  that  1^=1  can  imply  several  things :  ( 1 )  a  harmonic 

bifurcation,  (2)  a  periodic  turning  point,  and  (3)  integration  for 
one  period.  It  is  sufficient  to  note  that  these  harmonic 
bifurcations  exist  and  that  they  produce  unstable  periodic  solutions 
(Figure  9b);  they  are  not  explored  further  here  since  these 
bifurcations  do  not  indicate  toral  solutions. 

Keeping  the  same  value  of  Rag  as  in  Figure  7,  and  varying  the 


magnitude  of  a,  we  can  see  how  the  locations  of  the  Hopf  bifurcation 
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points  and  Hj,  the  harmonic  bifurcation  points  G^,  the  periodic 
turning  points  Gks-  and  the  diffeo-Hopf  points  Dj  and  Dffi  change  on  Pffl 
and  Pj  (Figure  10).  We  find  that  the  (large  dots  on  curve)  and  Hj 

(solid  curve)  curves  intersect  to  create  at  a^  =  0.49343,  as 

calculated  from  (2.49).  Parallel  to  these  Hopf  curves  are  the  curves 
for  the  periodic  turning  points  (small  dot)  and  G^s  (large  dash) 

that  represent  the  ends  of  the  periodic  branches.  These  curves  also 
intersect  at  a^j ,  indicating  that  the  periodic  branches  end  at  the 

same  value  of  Ra^..  As  the  branches  are  separated,  the  lines  also 

separate;  the  end  of  the  j  branch  becomes  farther  from  that  of  the  m 
branch  when  a  <  a^j .  The  reverse  is  true  for  a  >  a^. 

From  Figure  10,  we  find  that  there  are  only  single  diffeo-Hopf 

lines  on  each  side  of  a  •:  D-  (small  dash)  on  P-  when  a  <  a-,  and 

mj  j  J  mj 

Dm  (dash  dot)  on  Pffl  when  a  >  a^.  In  both  cases,  the  diffeo-Hopf 

points  occur  on  the  periodic  solution  that  branches  second  from  the 
conductive  solution.  For  a  stable  toral  solution,  as  occurs  in 
rotating  convection  (Wells  and  Shirer,  1987),  a  sufficient  condition 
is  that  both  diffeo-Hopf  curves  occur  to  the  right  of  with 

before  Dj,  or  the  left  of  a^j  with  Dj  before  Dm.  With  only  a  single 

curve  of  diffeo-Hopf  points  on  one  side  of  for  a  given  we 


Aspect  Ratio 


Figure  10  Variation  of  the  values  of  the  bifurcation  points  with 
aspect  ratio  when  Rag= 1.6836,  m=2,  j=l ,  P=7  and  r=0.01. 

The  periodic  solutions  grow  from  the  conductive 
solution  at  (curve  through  large  dots)  and 

(solid)  and  end  at  the  periodic  regular  turning  points 
(small  dots)  and  Gjg  (large  dash).  The  unstable 

toral  region  (shaded)  begins  at  the  diffeo-Hopf  points 
Dm  (dash  dot)  and  Dj  (small  dash)  and  extends  toward 

larger  values  of  Ra^.;  the  actual  edge  of  this  region  is 
not  known. 
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see  that  the  toral  solution  is  unstable.  If  the  diffeo-Hopf  points 
were  to  occur  on  the  first  branching  periodic  solution,  then  the 
toral  solution  would  be  stable.  This  would  occur  if  the  curve  for  D 

m 

sloped  to  the  left  rather  than  to  the  right  of  a^.  Thus,  by  knowing 
the  slope  of  these  diffeo-Hopf  curves  at  a  =  a^  and  then  by 

comparing  them  to  the  slopes  for  the  Hopf  curves,  we  can  infer  the 
stability  information  of  the  branching  quas i -periodic  solutions  in 
the  entire  toral  region. 

At  Rag  =  5.0  (Figure  11),  the  1^,  Hj,  G^,  and  Dj  curves  are 

similar  to  those  for  the  case  Rag  =  1.6836  in  Figure  10.  Now,  there 

are  two  additional  curves  representing  the  harmonic  bifurcations  Gme 

(small  dots  on  curve)  and  G’„  (long  dash),  at  which  the  periodic 

branches  lose  stability.  As  the  magnitude  of  Rag  is  increased, 

changes  in  the  size  of  the  toral  region  depend  on  how  far  Dj  is  from 

Hj,  arai  also  on  how  far  Dffl  is  from  HJn,  or  equivalently  on  the 

distance  D  •  and  D  have  moved  up  the  periodic  branches .  In  the  range 
j  m 

of  Rag  -  1.6836  to 
remains  about  constant,  since  the  value  of  |i^1l  at  Dj  In  this  range 

of  Rag  is  about  constant  at  =  1*2  (Figures  7,  8  and  12).  The 


too  -  65.0,  the  size  of  the  unstable  toral  region 


size  of  the  toral  region  also  depends  of  the  range  of  a,  since  at 
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Figure  11  Variation  of  the  values  of  the  bifurcation  points  with 

aspect  ratio  when  Ra^=5 . 0 ,  m=2,  J=1 ,  P=7,  and  rrO.Ol.  The 

periodic  solutions  PB  and  Pj  grow  from  the  conductive 
solution  at  H,,,  (curve  through  large  dots)  and  Hj  (solid), 

lose  stability,  or  gain  a  positive  characteristic 
multiplier  if  already  unstable,  at  (curve  through 

small  dots)  and  Gjg  (long,  thin  dash),  and  end  at  the 
periodic  regular  turning  points  Gm-  (small  dots)  and  Gj- 

( large  dash).  The  unstable  toral  region  (shaded)  begins 
at  the  diffeo-Hopf  points  Dm  (dash  dot)  and  Dj  (small 

dash)  and  extends  to  larger  values  of  Ray;  the  actual  edge 
of  this  region  is  not  known. 
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larger  values  of  Rag,  the  Dm  and  Dj  curves  extend  to  larger  values  of 

a.  The  regions  of  unstable  tori  are  shaded  in  Figures  10  and  11. 

For  the  values  of  Rag  i  65.0  that  have  been  investigated,  we  have  not 

found  any  stable  tori,  apparently  implying  that  stable  toral 
solutions  do  not  characterize  double  diffusive  convection.  We  can 
conclude  that  the  existence  of  a  double  Hopf  point  does  not  always 
guarantee  the  existence  of  stadole  tori  for  nearby  values  of  a,  but  it 
does  imply  that  there  will  be  either  stable  or  unstable  tori. 

Although  the  values  of  Ra^,  and  Rag  used  in  Figures  10  and  11 

give  steady  and  periodic  solutions,  the  values  of  the  forcing  rates 
are  far  from  the  realistic  ones  represented  by  the  Rag  =  Retp  line  of 

Figure  3.  Using  these  aphysical  values  Rag  =  1.6836  and  Ratp  =  9.7  of 

the  forcing  rates  produces  an  unrealistic  background  density  profile. 
This  is  evident  in  the  steady  solution  shown  in  Figure  13. 

Since  toral  solutions  involve  the  S„ _  and  T  terms  of  the 

modifications  of  the  background  salinity  and  temperature  profiles, 
and  since  the  toral  solutions  are  unstable,  these  terms  are  not 
activated.  Thus  the  modifications  of  background  density  profiles  cam 
be  expressed  exclusively  in  terms  of  the  TQ2  and  SQ2  coefficients. 

Although  the  background  temperature  and  salinity  modifications 

S02sin(2z)  and  T02sin(2z)  from  the  Fourier  expansions  of  S*  and  T* 
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act  to  decrease  the  interior  gradients  as  expected,  the  background 
density  profile  does  not  become  constant  with  height  because  the 
addition  of  pQ2 sin(2z)  to  pQQ *■  pq( z)  is  negligible  (Figure  13c). 

Because  of  the  fact  that  Ra^  >>  Rag,  which  corresponds  to  a  large 

amount  of  heating  from  below  destabilizing  the  system,  the  background 
density  increases  rapidly  with  height.  This  shows  that  straying  only 
somewhat  from  the  constant  density  line  (Figure  3)  can  lead  to 
aphysical  solutions. 

After  being  cautioned  as  to  the  care  required  in  choosing  values 
of  Rag  and  Ra^,,  we  would  like  to  know  physically  why  the  Pm  and  Pj 

solutions  exist,  why  the  steady  solution  is  preferred,  and  why  the 
periodic  solutions  end  abruptly.  These  topics  are  discussed  in  the 
next  section. 

4.2  Physical  Interpretation  of  the  Periodic  and  Steady  Solutions 

We  begin  by  finding  values  of  the  parameters  in  the  definitions 
(2.30)  and  (2.31)  of  Rag  and  R&p  that  satisfy  the  constant  density 

relation  PAS  =  <*AT.  The  density  ratio  aAT/PAS  represents  a  very 
important  observable  parameter  in  the  ocean.  For  example,  large 
values  of  oAT/PAS  imply  r  greater  destabilization  by  heat  than 
stabilization  by  salt.  This  density  ratio  will  prove  useful  in 
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guiding  our  examination  of  thermohaline  convection  in  estuaries  and 
under  melting  sea  ice. 

4.2.1  Thermohaline  Convection  in  Estuaries 

Estuaries  are  areas  of  interaction  between  fresh  and  salty 
water,  where  a  river  current  meets  the  ocean.  Here,  fresh  water 
overrides  salty  water,  which  yields  the  desired  salinity  gradient  for 
thermohaline  convection.  The  large  difference  in  salinity  between 
the  fresh  and  salty  water  produces  a  density  gradient  that  drives  a 
circulation.  Although  the  density  of  the  water  depends  on  both 
temperature  and  salinity,  in  the  estuary  case ,  the  temperature  range 
is  generally  small  (Dyer,  1973). 

As  an  example,  the  longitudinal  distributions  of  temperature  and 
salinity  in  the  Hardangerf jord  in  Norway  are  shown  in  Figure  14.  In 
the  winter,  cool,  fresh  water  overlies  warm,  salty  water.  Setting 

the  salinity  and  temperature  at  the  top  equal  to  33%.  and  5.5°C,  and 
at  the  bottom  to  35%.  and  6.5°C  gives  AS=2%.  and  hT=l°C.  Also,  we 
set  TQ0=  280  K,  and  S0Q=35%..  These  values  give  aAT//?AS  =  6.25xl0'a, 

which  is  far  from  one  and  the  constant  density  line.  Also,  these 
parameter  values  give  Rag  -  4.024  and  Ra^.  =  0.252.  However,  these 

values  of  the  forcing  rates  are  well  below  those  for  the  constant 
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density  and  Hopf  lines  in  Figure  3,  thereby  indicating  that  a 
conductive  solution  would  be  stable.  The  salt  gradient  is  so 
stabilizing  that  convection  cannot  occur.  We  conclude  that  the 
wintertime  estuary  case  does  not  yield  thermohaline  convection. 

4.2.2  Thermohaline  Convection  Under  Melting  Sea  Ice 

Melting  sea  ice  provides  a  flux  of  fresh  water  to  the  surface 
layer  of  the  ocean,  possibly  producing  the  desired  salinity  gradient 
for  thermohaline  convection.  In  order  to  have  periodic  solutions  and 
a  realistic  background  density,  we  choose  the  values  of  Rag  near  that 

for  RaT  and  to  the  right  of  the  intersection  of  the  two  curves  in 
Figure  3.  We  will  choose  Rag  =  65.0  and  Rap  =  64.94.  The 
bifurcation  diagram  for  Rag  =  65.0  is  shown  in  Figure  12,  which  is 

similar  to  the  cases  in  Figures  7  and  8. 

Upon  knowing  the  values  of  the  forcing  rates,  we  find  that 
choosing  the  parameter  values  for  the  sea  ice  case  is  quite 
straightforward.  We  choose  TQ0=  283  K,  the  temperature  of  the 

sumner time  polar  seasonal  thennocline  (Figure  15c).  The  temperature 
at  the  top  must  be  near  that  of  melting  ice,  277  K.  Choosing  a  lower 
temperature  is  inconsistent  with  the  linear  equation  of  state  (2.14) 
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(«)  (b)  (c) 


Tropics  Mid  latitudes  Polar 


Figure  15  Schematic  temperature  distributions  with  depth  for 

tropical  (a),  mid-latitude  (b),  and  polar  regions  (c) 
for  simmer  (solid)  and  winter  (dashed)  conditions.  The 
major  regimes  in  the  vertical  are  the  near-surface 
mixed  layer,  the  seasonal  thennocline,  and  the  deep, 
permanent  thennocline  (from  Apel,  1987). 
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used  to  develop  the  model.  Thus ,  we  obtain  AT  =  6  K.  As  before,  we 
specify  SQ0=  35%  «. 

In  order  to  find  the  values  of  the  other  parameters  and  the 
magnitudes  of  the  forcing  rates  from  (2.30)  and  (2.31),  we  choose  a 
value  for  v.  We  note  that  for  large-scale  phenomena,  the  horizontal 
eddy  diffusion  term  dominates  (Table  2).  However,  on  the  convective 

scale  i/y  -  t'k  ~  100  ma/s,  which  is  the  smallest  value  in  the  range  of 

given  in  Table  2.  Since  P  =  I'/kp,  we  set  =  14.286  m2/s  in 

order  to  have  P  =  7. 

Next,  we  must  choose  AS  and  Zj  so  that  ka^  =  Raj.  The  value  of 

AS  for  which  AT/Too  =  AS/S0Q  is  AS  =  0.742%..  When  Ra^  =  64.94,  we 
find  from  the  value  of  v  and  (2.31)  that  Zj  ~  351.7  m.  This  depth 
agrees  relatively  well  with  that  of  the  s miner  polar  seasonal 


Table  2  Typical  values  of  horizontal  ard  vertical  eddy  diffusion 
coefficients  for  large-scale  motion  of  sea  water  (from 
Apel,  1987). 


=  3xl0~5  to  2x10  2  m*/s 
vh  =  10®  to  10s  ma/s 
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thennocline  (Figure  15).  Also,  this  depth  is  verified  in  Figure  16 
by  the  location  of  the  asterisk.  Here,  we  combine  (2.48)  with  (2.30) 
and  (2.31)  and  let  a 
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This  equation  gives  values  of  Zj  that  are  near  the  constant  density 
line  and  the  double  Hopf  point. 

The  spectral  components  for  the  periodic  solution  Pm  for  Rag  = 


65.0  and  Ra^  =  64.94  are  shown  in  Figure  17.  We  note  that  the 

salinity  field  (coefficients  Smi  and  SQ2 )  lags  the  temperature  field 

(coefficients  T  and  T  _ )  because  of  the  differences  in 
mi  02 

diffusivities  represented  by  the  ratio  t,  as  discussed  in  Chapter  2. 
However,  to  understand  what  the  interaction  between  these  components 
means  physically,  we  must  convert  them  to  their  dimensional  forms 
using  (2.21)— (2. 26 ) .  Although  the  model  gives  numerical  values  for 
the  temperature  and  salinity  perturbations,  it  is  difficult  to 

9  9 

explain  the  characteristics  of  the  circulation  unless  T  and  S  are 
combined  to  yield  the  perturbation  density.  Combining  (2.15),  (2.7), 
and  (2.8),  we  f ind  that 
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Figure  16  Values  of  domain  depth  (m)  for  the  double  Hopf  point 

with  the  following  values:  1^14.286  a*/s,  v=100  aa/s, 
P=7,  Toq=283  K,  S00=35*.,  AT=6  K,  AS =0.742*.,  m=2,  j=l, 

g=9.8  s/s* ,  and  8^=0.4934.  Ihe  constant  density  line 

is  dashed,  and  the  asterisk  denotes  the  sea  ice  case 
discussed  in  the  text. 


and 
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Substituting  Sb  =  SQ2sin(2z*>  and  =  TQ2sin(2z*),  where  b  indicates 


the  background  state ,  into  (2.16)  gives 
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where  SQ2  and  TQ2  are  calculated  from  solutions  to  the  model. 
Substituting  Sc  =  Smicos(mx*)sin(z*)  and  Tc  =  Tmicos(mx* )sin( z* ) , 


where  c  donotes  the  convective  state,  into  (2.16)  gives 


p  S  T 

in  -  _aii  _  iffii 

poo  Soo  Too 


(4.8) 


where  S_  and  T  .  are  also  calculated  from  the  solutions, 
nu  mi 

The  variations  with  time  of  the  variables  Pmi  and  pQ2,  along 

!  r 

with  y  ,  are  shown  in  Figure  18.  Since  p  w  >0,  where  the  overbar 

denotes  a  horizontal  average,  implies  an  indirect  circulation,  for 
which  more  dense  fluid  is  transported  up  and  less  dense  fluid  is 

I  ♦  ♦ 

brought  down,  we  let  pc  s  pmicos(mx  )sin(z  )  and 

w’  =  nr*  cosdnx*  >sin(  z* )  to  find  that  P-.Vl,.  >  0  also  implies  an 
mi  mi  mi 

indirect  circulation.  From  Figure  18,  we  see  that  the  flow  is 
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Figure  18  Variations  in  density  coefficients  pQZ  and  Pmi  and 

stream  function  coefficient  corresponding  to  the 

periodic  solution  in  Figure  17.  The  shaded  area 
represents  times  in  which  the  circulation  is  indirect. 
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characterized  briefly  by  an  indirect  circulation  between 
approximately  t  =  416.04  and  t  =  416.08,  and  again  between 
approximately  t  =  416.33  and  t  =  416.37.  This  indirect  circulation 
occurs  for  much  less  than  half  of  the  period.  This  behavior  is 
similar  to  that  seen  in  rotating  convection  (Pyle,  1986). 

To  examine  the  periodic  solution  in  more  detail,  we  contour  the 

t-  1  $  < 

values  of  rc  =  fmisin(mx  )sin(z  )  and  pc  =  Pmicos(mx  )sin(z  )  as 

functions  of  x  and  z  ,  and  we  plot  p^  =  -Ap  z/Zj,  +  pQ2 sin(2z  ) 

versus  height  (Figures  19a- 19g)  at  various  times  throughout  a  single 
period  (Figures  17  and  18).  At  t  =  416.00  (Figure  19a),  the 
circulation  is  direct,  given  by  Pmivmi  <  0:  there  is  an  upward 

transport  of  lighter  fluid  and  a  downward  transport  of  heavier  fluid. 
The  background  density  increases  with  height  implying  static 

instability  in  the  region  of  0  <  z*  <  w/4  and  3n/4  <  z*  <  ",  and  the 
fluid  responds  by  attempting  to  destroy  these  unstable  gradients. 

The  contribution  to  the  total  background  density  by  Ap  in  pQ  is 

negligible,  since  this  case  is  almost  on  the  constant  density  line. 

Although  less  dense  fluid  is  moving  up,  and  there  is  a  direct 
circulation,  the  stratification  is  unstable  near  the  boundaries.  The 
only  way  to  eliminate  the  unstable  density  profiles  at  the  boundaries 
is  by  having  an  indirect  circulation  (Pmi^mi  >  0)  that  moves  more 
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Figure  19  (cont. )  (e)  At  tine  t=416.36,  the  circulation  is  indirect  since 

p  has  become  positive.  The  background  density  in 
nil 

unstable  in  the  interior  as  is  now  negative. 
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dense  fluid  upwards  (Figure  19b).  This  indirect  circulation  is 
short-lived  since  it  creates  another  unstable  gradient  in  the  middle 
of  the  density  profile  that  lasts  longer  than  the  indirect 
circulation.  Note  that  the  indirect  circulation  always  occurs  just 
after  the  value  of  pqz  becomes  negative  (Figure  18),  implying  that  an 

unstable  interior  density  gradient  occurs.  Now,  a  direct  circulation 
attempts  to  destroy  this  interior  instability  by  bringing  upwards 
lighter  fluid  that  is  created  by  the  forcing  from  below  (Figure  19c). 
This  process  is  repeated  in  Figures  19e-g  to  complete  the  period. 

The  search  for  an  equilibrium  having  no  density  variation  with 
height  produces  a  periodic  solution,  since  an  attempt  to  destroy  an 
unstable  gradient  in  one  regime  produces  an  unstable  gradient  in 
another.  Eventually,  this  attempt  to  gain  equilibrium  breaks  down  as 
the  amplitudes  of  the  oscillations  in  the  background  modification 
become  too  large,  and  the  flow  can  no  longer  create  the  necessary 
indirect  circulations  to  maintain  a  periodic  state.  This  breakdown 
occurs  when  the  periodic  branch  ends,  and  the  steady  solution  becomes 
the  preferred  one.  The  P-  solution,  which  becomes  stable  via  the 

J 

toral  solution,  has  characteristics  similar  to  those  for  the  Pm 
solution. 

The  above  behavior  of  the  periodic  solutions  may  be  understood 
also  by  comparing  the  bifurcation  diagrams  in  Figures  7,  8  and  12. 
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Here,  the  periodic  solutions  occur  only  at  small  values  of  IV^I  »  so 

these  solutions  weakly  transport  heat  and  momentum  vertically.  As 
the  value  of  the  forcing  rate  is  increased,  greater  vertical  heat 
transport  is  required,  which  the  oscillatory  state  cannot  maintain, 
and  so  the  periodic  solution  ends.  The  lower  branch  of  the  steady 
solution  is  probably  unstable  for  the  same  reason.  The  stable  steady 
solution  provides  the  necessary  large  vertical  transport  of  heat  to 
maintain  an  equilibrium.  Veronis  (1968)  found  that  for  P  i  1,  the 
oscillatory  solutions  were  replaced  with  steady  solutions,  that 
had  larger  heat  fluxes. 

An  example  of  a  steady  solution  S  for  the  same  Rag  as  in  the 

periodic  solution  discussed  above  is  shown  in  Figure  20.  As  before, 

) 

the  slope  owing  to  Ap  is  negligible,  and  the  contribution  from  p^ 

dominates  in  the  background  density  profile.  Here,  the  circulation 
is  always  a  direct  one  that  is  trying  to  alter  the  background  density 
gradient  in  the  interior,  where  more  dense  fluid  overlies  less  dense 
fluid,  in  order  to  create  a  background  density  profile  that  is 
constant  with  height.  The  steady  solution  does  not  accomplish  this 
task  because  the  forcing  rates  Rag  >  0  and  Ra^  >  0  produce  a  steady 

source  of  fresh  water  above  and  warm  water  below.  This  type  of 


forcing  produces  less  dense  fluid  at  the  boundaries,  which  can  be 
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Figure  20  Contours  of  convective  density  pQ  (kg/m  )  and  stream  function 
i>c  (m*/s),  and  the  modification  of  the  background  density  ^ 
(kg/m3 )  for  the  steady  solution  S#  with  the  values  Rag=65 , 
88^65.1,  K^H.286  m*/s,  v=100  ma/s,  P=7,  T0#=283  K, 

S00r35X.f  AT=6  K,  ASrO.742%.,  m=2,  j=l,  g=9.8  m/s*f  8^=0.493 

and  rsO.Ol.  Here,  the  circulation  is  always  direct,  and  the 
background  density  is  instable  in  the  interior. 
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seen  in  the  background  density  profile  (Figui'e  20). 

From  the  periodic  and  steady  solutions  in  Figures  19  and  20,  we 
see  that  if  the  forcing  rates  are  chosen  according  to  Figure  3,  hen 
realistic  density  perturbations  and  circulations  result.  Although  no 
stable  toral  solutions  are  found  in  the  region  of  physically 
realistic  solutions  owing  to  the  relationships  between  the  slopes  of 
the  diffeo-Hopf  and  Hopf  lines,  they  may  exist  for  some  ap».ysicai 
values  of  Rag  and  Ra^  not  examined  here.  The  above  results  and  their 

implications  are  summarized  in  the  concluding  chapter. 
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Chapter  5 
CONCLUSION 

Density  profile  modifications  caused  by  thermohal ine  convection, 
which  might  be  observed  in  estuaries  and  near  melting  sea  ice,  may 
lead  to  the  trapping  of  acoustical  energy  through  the  development  of 
a  near-surface  sound  channel  (Apel,  1987).  This  channel  exists  where 
the  density  decreases  with  height,  and  its  boundaries  are  defined  by 
unstable  density  gradients  in  the  same  way  temperature  inversions  act 
as  boundaries  in  the  atmosphere.  In  this  channel,  the  competition 
between  the  stabilizing  salinity  gradient  and  the  destabilizing 
temperature  gradient  allows  steady  and  periodic  solutions,  and  it  is 
these  solutions  that  modify  the  background  density  profiles.  We  have 
studied  these  modifications  by  examining  the  existence  of  Hopf  points 
and  their  associated  periodic  solutions.  We  have  also  studied  the 
existence  of  double  Hopf  points  and  diffeo-Hopf  points  as  well  as 
their  associated  toral  solutions.  Since  double  Hopf  points  exist  in 
the  partial  differential  equations  representing  thennohaline 
convection,  we  would  expect  that  toral  solutions  might  cause 
significant  modifications  of  the  background  density  profiles. 

However,  using  a  sixteen-coefficient,  nonlinear  spectral  model,  we 
found  that  these  toral  solutions  are  unstable  and  so  do  not  lead  to 
observable  changes  in  these  profiles.  Although  these  unstable  toral 
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solutions  do  not  modify  the  imposed  density  gradients,  they  do  cause 
the  otherwise  unstable  periodic  solution  of  wavenumber  j  to  gain 
stability.  We  conclude  that  the  presence  of  a  double  Hopf  point  does 
not  always  imply  the  existence  of  stable  toral  solutions,  while  the 
presence  of  a  diffeo-Hopf  point  does  always  imply  the  existence  of  a 
(stable  or  unstable)  toral  solution. 

In  contrast  to  the  values  used  by  Veronis  (1968),  Knobloch  and 
Proctor  (1981),  Knobloch  et  al.  (1986),  and  Huppert  and  Moore  (1976) 
to  study  thermohaline  convection  —  Prandtl  number  P  =  1  and  ratio  of 
diffusivities  r  =  0.125  to  0.5  —  we  use  the  parameter  values 
representing  sea  water  —  P  =  7  and  r  =  0.01.  Upon  estimating  the 
parameter  values  characterizing  circulations  in  estuaries,  we  find 
that  the  values  of  the  forcing  rates  admit  only  the  conductive 
solution.  This  result  occurs  because  of  the  large  difference  in  the 
values  of  the  solutai  and  thermal  forcing  rates:  Rag  >>  Rap. 

However,  because  Rag  =  Rap  in  the  melting  sea  ice  case,  the  rates  of 

forcing  produce  plausible  periodic  and  steady  solutions  that 
significantly  modify  the  background  density. 

In  the  periodic  regime  Rag  =  Rap  of  thermohal ine  convection,  the 
fluid  is  characterized  by  a  direct  circulation  for  a  majority  of  the 


period,  during  which  there  is  an  upward  transport  of  lighter  fluid. 
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The  fluid  is  trying  to  eliminate  the  unstable  density  gradient  that 
is  caused  by  the  forcing  near  the  boundaries.  By  attempting  to 
destroy  these  unstable  gradients,  the  flow  creates  another  unstable 
gradient  in  the  interior  of  the  fluid.  This  gradient  is  produced  by 
an  indirect  circulation  that  moves  denser  fluid  upward  and  lighter 
fluid  downward.  Therefore,  in  a  periodic  solution,  an  attempt  to 
destroy  an  undesirable  density  gradient  creates  another.  In  the 
periodic  case,  the  size  of  the  sound  channel  alternately  decreases 
and  increases  as  the  unstable  density  gradient,  which  does  not  allow 
acoustic  propagation,  alternately  moves  toward  and  away  from  the 
boundaries.  When  the  interior  density  gradient  is  stable,  there  is 
an  interior  sound  channel.  As  this  interior  gradient  becomes 
unstable,  the  sound  channel  separates  into  two  pieces  associated  with 
the  stable  gradients  at  the  boundaries.  Thus,  the  circulation  causes 
the  location  of  the  sound  channel  to  fluctuate  significantly  when  the 

background  density  is  near  neutral  so  that  Rag  =  Ra^- 

As  the  forcing  rate  is  increased,  the  magnitude  of  the 
background  density  modification  increases,  and  so  portions  of  the 
density  profiles  become  increasingly  more  unstable.  The  periodic 
solution  eventually  breaks  down  since  it  can  no  longer  transport 
fluid  indirectly  during  an  oscillation.  In  this  case,  the  periodic 
solutions  end  abruptly  through  periodic  regular  turning  points  for 
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all  values  of  the  solutal  forcing  rate  Rag  that  were  investigated. 

The  preferred  solution  is  then  a  steady  one,  which  is  a  direct 
circulation.  This  circulation  supports  an  unstable  gradient  in  the 
interior  of  the  fluid.  Thus,  the  sound  channel  is  found  in  two 
pieces  in  the  steady  case,  since  its  boundaries  are  now  the  top  and 
bottom  of  the  fluid  and  those  for  the  interior  unstable  gradient.  In 
the  steady  case,  there  is  always  a  near-surface  sound  channel. 

To  study  the  depth  and  location  of  the  sound  channel,  we  have 
investigated  periodic  and  steady  solutions.  The  steady  solutions  of 
horizontal  wavenumber  m  and  j  are  found  analytically,  and  the 
horizontal  wavenumber  m  and  j  periodic  solutions,  Pm  and  Pj  are  found 

numerically,  as  the  thermal  forcing  rate  RaT  is  increased  past  that 
for  the  wavenumber  m  and  j  Hopf  bifurcation  points,  and  Hj.  At  a 
=  amj ,  and  Hj  coincide  to  form  the  double  Hopf  point  .  For  an 
aspect  ratio  a  <  a^j,  separates  into  and  Hj.  and  the  Pm 
solution  branches  first;  for  a  >  a^j,  Pj  branches  first.  An 
eigenvalue  analysis  of  the  Pffl  branch,  when  a  <  a^,  reveals  that  a 
harmonic  bifurcation  occurs  for  some  values  of  Rag,  but  that  no 

diffeo-Hopf  points  exist  from  which  toral  solutions  would  emanate 
from  the  periodic  solution.  For  values  just  past  that  for  H ■ ,  the  P^ 

J  v 


95 


branch  is  unstable,  but  gains  stability  through  a  diffeo-Hopf 
bifurcation  It  is  this  bifurcation  that  leads  to  an  unstable 

torus.  With  the  existence  of  only  one  diffeo-Hopf  point  on  the 
unstable  periodic  solution,  there  is  no  stable  secondary  connecting 
branch  linking  the  wavenumber  m  and  j  periodic  branches. 

Because  of  the  unstable  toral  solution,  there  is  no  exchange  of 
stability  between  the  horizontal  wavenumber  m  and  j  periodic 
solutions,  as  was  found  for  steady  convection  by  Chang  and  Shirer 
(1984)  and  for  rotating  convection  by  Wells  and  Shirer  (1987).  All 
energy  put  into  the  catalytic,  or  wavenumber  interaction  terms 
(m+j , 2 )  and  (m-j,2)  dissipates  quickly.  Also,  all  energy  in  the 
terms  having  horizontal  wavenumber  zero  and  vertical  wavenumber  four 
that  represent  the  modification  of  the  background  density  profile  by 
the  toral  solution  damps  to  zero.  Thus,  we  are  able  to  represent  the 
modification  of  the  background  density  field  by  the  (0,2)  terms 
alone . 

As  the  value  of  a  is  varied  from  that  for  a^,  the  values  of  the 

diffeo-Hopf  points  Ra^(Dj)  and  Rarp(Dm)  on  the  periodic  solutions  Pj 

and  P  also  vary,  and  these  curves  form  the  boundary  of  the  unstable 
m 

toral  region.  With  only  a  D;  bifurcation  point  on  the  unstable  P^ 

J  J 

solution  for  a  <  a^  and  a  Dm  point  on  the  unstable  solution  Pffi  for  a 


>  ,  we  have  only  a  single  diffeo-Hopf  line  on  each  side  of  . 

The  diffeo-Hopf  point  is  always  on  the  secorxi  periodic  branch,  so 
that  R&plHjjj)  <  Rap(Hj  )  <  Rap(Dj  )  for  a  <  ,  and  R&p(Hj )  <  Rap(Hm)  < 

Rap(Dm)  for  a  >  a^j ,  and  so  there  is  always  an  unstable  toral 

solution.  Therefore,  by  knowing  the  relationships  between  the  Hopf 
and  diffeo-Hopf  lines,  we  can  determine  whether  the  bifurcating 
quasi-periodic  solution  is  stable  or  unstable. 

In  rotating  convection,  where  Dm  and  both  existed  on  one  side 

of  the  double  Hopf  point  so  that  Raq,(Hm)  <  RaT(Dm)  <  Ra^,(Hj)  < 

RaT(Dj)  (Wells  and  Shirer,  1987),  stable  toral  solutions  were  found. 

The  truncation  chosen  in  the  present  sixteen-coefficient  model  is 
similar  to  that  used  by  both  Chang  and  Shirer  (1984)  and  Wells  and 
Shirer  ( 1987 ) ,  whose  models  successfully  captured  the  observed 
transfer  of  energy  from  larger  to  smaller  horizontal  wavenunbers . 

The  relationships  between  the  diffeo-Hopf  curves,  which  bovnd 
the  toral  region,  and  the  Hopf  curves  can  be  approximated  using  a 
power  series  analysis  of  the  partial  differential  equations.  This 
analysis,  which  is  still  being  developed,  will  also  provide  the 
truncation  level  required  for  modeling  toral  solutions  and  the 
relative  importance  of  each  component  of  the  truncation  in  the 


representation  of  these  toral  solutions.  With  this  analysis,  we  will 
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be  able  to  see  whether  the  trigonometric  representation  of  the 
variables  in  the  present  model  is  sufficient,  and  how  sensitive  the 
existence  and  stability  of  the  toral  solution  are  to  the  nimher  of 
terms  included  in  the  expansion. 
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